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

/*
 * wmmarr.cpp
 * arrays of WMM modes, mode interference evaluation
 */

#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"wmmvis.h"
#include"mode.h"
#include"wmmfld.h"
#include"wmmarr.h"

/* WMM_ModeArray: Array of WMM_Modes */

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

/* initialize */
WMM_ModeArray::WMM_ModeArray()
              : num(0)
{
	int i;
	for(i=0; i<=MAXNUMMODES-1; ++i) mvec[i] = NULL;
}

/* destroy */
WMM_ModeArray::~WMM_ModeArray()
{
	int i;
	for(i=0; i<=MAXNUMMODES-1; ++i) 
	{
		if(mvec[i] != NULL)
		{
			delete mvec[i];
			mvec[i] = NULL;
		}
	}
	num = 0; 
}

/* copy constructor */
WMM_ModeArray::WMM_ModeArray(const WMM_ModeArray& ma)
	      : num(ma.num)
{
	int i;
	for(i=0; i<=num-1; ++i)
	{
		mvec[i] = new WMM_Mode;
		*mvec[i] = *(ma.mvec[i]);
	}
}
	
/* assignment */
WMM_ModeArray& WMM_ModeArray::operator=(const WMM_ModeArray& ma)
{
	int i;
	if(this != &ma)
	{
		for(i=0; i<=MAXNUMMODES-1; ++i) 
		{
			if(mvec[i] != NULL)
			{
				delete mvec[i];
				mvec[i] = NULL;
			}
		}
		num = ma.num; 
		for(i=0; i<=num-1; ++i)
		{
			mvec[i] = new WMM_Mode;
			*mvec[i] = *(ma.mvec[i]);
		}
	}
	return *this;
}

/* delete all Mode entries */
void WMM_ModeArray::clear()
{
	int i;
	for(i=0; i<=MAXNUMMODES-1; ++i) 
	{
		if(mvec[i] != NULL)
		{
			delete mvec[i];
			mvec[i] = NULL;
		}
	}
	num = 0; 
}

/* input from FILE dat */
void WMM_ModeArray::read(FILE *dat)
{
	int i;
	if(dat == stderr || dat == stdout)
	{
		wmmarrerror("WMM_ModeArray: read");
	}
	else
	{
		for(i=0; i<=MAXNUMMODES-1; ++i) 
		{
			if(mvec[i] != NULL)
			{
				delete mvec[i];
				mvec[i] = NULL;
			}
		}
		num = 0; 
		num = fgetint(dat);
		if(num >= MAXNUMMODES)
		{
			fprintf(stderr, 
				"WMM_ModeArray: read: NUM reduced to %d.\n",
				MAXNUMMODES);
			num = MAXNUMMODES;
		}
		for(i=0; i<=num-1; ++i)
		{
			mvec[i] = new WMM_Mode;
			(*mvec[i]).read(dat); 
		}
	}
	return;
}

/* output to FILE dat */
void WMM_ModeArray::write(FILE *dat)
{
	int i;
	comment("WMM_ModeArray", dat);
	comment("num", dat); fputint(num, dat);
	for(i=0; i<=num-1; ++i)
	{
		comment("*** -------- ***", dat);
		(*mvec[i]).write(dat); 
	}
	return;
}

/* input from default file */
void WMM_ModeArray::read_def(char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	int i;
	name[0] = 'w';
	name[1] = 'm';
	name[2] = 'm';
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 'm';
	name[7] = 'a';
	name[8] = 'r';
	name[9] = 0;

	fprintf(stderr, "<< %s ", name);
	
	dat = fopen(name, "r");
	read(dat);
	fclose(dat);

	fprintf(stderr, "#%d\n", num);
	for(i=0; i<=num-1; ++i)
		fprintf(stderr, 
		  "   (%d) %c%c%c,  beta: %.10g,  neff: %.10g,  npcB: %.10g\n",
			i,
			poltochar1((*mvec[i]).pol),
			poltochar2((*mvec[i]).pol),
			symtochar((*mvec[i]).sym),
			(*mvec[i]).beta, 
			(*mvec[i]).neff, 
			(*mvec[i]).npcB); 
	return;
}

/* output to default file */
void WMM_ModeArray::write_def(char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	int i;
	name[0] = 'w';
	name[1] = 'm';
	name[2] = 'm';
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 'm';
	name[7] = 'a';
	name[8] = 'r';
	name[9] = 0;

	fprintf(stderr, ">> %s #%d\n", name, num);
	for(i=0; i<=num-1; ++i)
		fprintf(stderr, 
		  "   (%d) %c%c%c,  beta: %.10g,  neff: %.10g,  npcB: %.10g\n",
			i,
			poltochar1((*mvec[i]).pol),
			poltochar2((*mvec[i]).pol),
			symtochar((*mvec[i]).sym),
			(*mvec[i]).beta, 
			(*mvec[i]).neff, 
			(*mvec[i]).npcB); 
	dat = fopen(name, "w+");
	write(dat);
	fclose(dat);

	return;
}

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

/* add a mode */
void WMM_ModeArray::add(WMM_Mode m)
{
	if(num >= MAXNUMMODES)
		wmmarrerror("WMM_ModeArray: add: MAXNUMMODES exceeded");
	if(mvec[num] != NULL)
		wmmarrerror("WMM_ModeArray: add: mvec corrupted");
	mvec[num] = new WMM_Mode;
	*mvec[num] = m;
	++num;
	return;
}

/* delete a mode entry */
void WMM_ModeArray::remove(int i)
{
	int j;
	if(i<0 || i>=num) 
		wmmarrerror("WMM_ModeArray: remove: argument out of range");
	if(mvec[i] == NULL)
		wmmarrerror("WMM_ModeArray: remove: mvec corrupted");
	delete mvec[i];
	mvec[i] = NULL;
	if(i==num-1) 
	{
		--num;
		return;
	}
	for(j=i; j<=num-2; ++j) mvec[j] = mvec[j+1];
	mvec[num-1] = NULL;
	--num;
	return;
}

/* add an entire WMM_ModeArray nma, nma is cleared ! */ 
void WMM_ModeArray::merge(WMM_ModeArray& ma)
{
	int j;
	if(num+ma.num >= MAXNUMMODES)
		wmmarrerror("WMM_ModeArray: merge: MAXNUMMODES exceeded");
	if(mvec[num] != NULL)
		wmmarrerror("WMM_ModeArray: merge: mvec corrupted");
	for(j=0; j<=ma.num-1; ++j)
	{
		if(mvec[num] != NULL)
			wmmarrerror("WMM_ModeArray: merge: mvec corrupted");
		mvec[num] = ma.mvec[j];
		++num;
		ma.mvec[j] = NULL;
	}
	ma.num = 0;
	return;
}

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

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

/* 
 * mode interference evaluation and visualization 
 * all modes are assumed to belong to the same waveguide ! 
 */

/* field superposition at point (x, y, z), 
   amp:  complex amplitudes at z=0
	 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, 
         cbet(j) = this(j).beta+pert(j)
   cp:   EX, EY, EZ, HX, HY, HZ, SZ */
Complex WMM_ModeArray::field(Cvector amp, Cvector pert, Fcomp cp, 
                             double x, double y, double z) const
{
	Complex s, iphase;
	double f, sum;
	Cvector a;
	int i, j, k;
	a = amp;	
	if(amp.nel < num)
	{
		a = Cvector(num);
		a.init(CC0);
		for(i=0; i<=amp.nel-1; ++i) a(i) = amp(i);
	}
	Cvector cbet(num);
	for(i=0; i<=num-1; ++i) cbet(i) = pert(i)+(*mvec[i]).beta;
	if(z != 0.0)
	{
		for(i=0; i<=num-1; ++i) 
		{
			iphase = CCI*cbet(i)*(-z);
			a(i) = a(i)*exp(iphase);
		}
	}
	s = CC0;
	if(cp != SZ)
	{
		for(i=0; i<=num-1; ++i)
		{
			f = (*mvec[i]).field(cp, x, y);  
			if(cp != EZ && cp != HZ) 
				s = s+a(i)*f;
			else
			{
				a(i) = CCI*a(i);
				s = s+a(i)*f;
			}
			
		}
		return s;
	}
	sum = 0.0;
	Dvector fEx(num); 
	Dvector fEy(num); 
	Dvector fHx(num); 
	Dvector fHy(num); 
	for(j=0; j<=num-1; ++j) fEx(j) = (*mvec[j]).field(EX, x, y);
	for(j=0; j<=num-1; ++j) fEy(j) = (*mvec[j]).field(EY, x, y);
	for(j=0; j<=num-1; ++j) fHx(j) = (*mvec[j]).field(HX, x, y);
	for(j=0; j<=num-1; ++j) fHy(j) = (*mvec[j]).field(HY, x, y);
	for(j=0; j<=num-1; ++j)
	{
		for(k=0; k<=num-1; ++k)
		{
			f = (a(j).re*a(k).re + a(j).im*a(k).im);
			sum += f*(fEx(j)*fHy(k)-fEy(j)*fHx(k));
		}
	}
	s = Complex(0.5*sum);
	return s;
}

/* evaluate component cp on a rectangular npx x npy mesh on area disp 
   at position z
   amp:  complex amplitudes at z=0
	 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, 
         cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ 
   foa: MOD, SQR, REP, IMP */
Dmatrix WMM_ModeArray::fieldmat(Cvector amp, Cvector pert, double z,
                                Rect disp, int npx, int npy, 
                                Fcomp cp, Afo foa) const
{
	double x, y;
	double dx, dy;
	int i, j;
	if(npx <= 1) wmmarrerror("fieldmat: npx <= 1");
	if(npy <= 1) wmmarrerror("fieldmat: npy <= 1");
	Dmatrix f(npx, npy);
	Complex ft;
	double ftd;
	dx = (disp.x1-disp.x0)/(npx-1);
	dy = (disp.y1-disp.y0)/(npy-1);
	for(i=0; i<=npx-1; ++i)
	{
		x = disp.x0+i*dx;	
		for(j=0; j<=npy-1; ++j)
		{
			y = disp.y0+j*dy;	
			ft = field(amp, pert, cp, x, y, z);
			ftd = 0.0; 
			switch(foa)
			{
			case MOD: case ORG: 
				ftd = ft.abs(); 
				break;
			case SQR: 
				ftd = ft.sqabs(); 
				break;
			case REP:  
				ftd = ft.re;
				break;
			case IMP:  
				ftd = ft.im;
				break;
			}
			f(i, j) = ftd;
		}
	}
	return f;
} 

/* store nump values of component cp between (x0, y0) and (x1, y1) in a vector 
   at position z
   amp:  complex amplitudes at z=0
	 amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, 
         cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ 
   foa: MOD, SQR, REP, IMP */
Dvector WMM_ModeArray::fieldvec(Cvector amp, Cvector pert, double z, 
				Fcomp cp, Afo foa, int nump, 
                                double x0, double y0, 
				double x1, double y1) const
{
	double x, y;
	double dx, dy;
	int j;
	if(nump <= 1) wmmarrerror("fieldvec: nump <= 1");
	Dvector f(nump);
	Complex ft; 
	double ftd;
	dx = (x1-x0)/(nump-1);
	dy = (y1-y0)/(nump-1);
	for(j=0; j<=nump-1; ++j)
	{
		x = x0+j*dx;	
		y = y0+j*dy;	
		ft = field(amp, pert, cp, x, y, z);
		ftd = 0.0;
		switch(foa)
		{
		case MOD: case ORG: 
			ftd = ft.abs(); 
			break;
		case SQR: 
			ftd = ft.sqabs(); 
			break;
		case REP:  
			ftd = ft.re;
			break;
		case IMP:  
			ftd = ft.im;
			break;
		}
		f(j) = ftd;
	}
	return f;
}

/* - Three segment coupler, power transfer evaluation -------------------- 
   imode: input mode, excites the modes this(j) with rel. power 1 at z=0 
   omode: output mode, relative power is returned 
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j) 
(!) this(j) are assumed to be normalized. */

/* weights for the power evaluation:
   w = ( omode | this(m) ) ( this(m) | imode ) */
double WMM_ModeArray::pweight(const WMM_Mode& imode, int m,
                              const WMM_Mode& omode) 
{
	return lscalprod(omode, (*mvec[m]))*lscalprod((*mvec[m]), imode);
}
   
/* single relative power level for device length l */    
double WMM_ModeArray::iopower(const WMM_Mode& imode, const WMM_Mode& omode, 
                              Cvector pert, double l)
{
	Complex a, am, cbet;
	int m;
	a=CC0;
	for(m=0; m<=num-1; ++m)
	{
		cbet = pert(m)+(*mvec[m]).beta;
		cbet = CCI*cbet*(-l);
		am = exp(cbet)*pweight(imode, m, omode);
		a = a+am;
	}
	return a.sqabs();
}

/* output power for numl devices of lengths between lmin and lmax */    
Dvector WMM_ModeArray::iopower(const WMM_Mode& imode, const WMM_Mode& omode, 
                               Cvector pert, 
                               int numl, double lbeg, double lend)
{
	Complex a, am, cbet;
	a=CC0;
	Dvector w(num);
	double l, dl;
	int li, m;
	if(numl <= 0) wmmarrerror("iopower: numl <= 0");
	Dvector p(numl);
	if(numl == 1) 
	{
		p(0) = iopower(imode, omode, pert, 0.5*(lbeg+lend));
		return p;
	}
	for(m=0; m<=num-1; ++m)
	{
		w(m) = pweight(imode, m, omode);	
	}
	dl = (lend-lbeg)/(numl-1.0);
	for(li=0; li<=numl-1; ++li)
	{
		a=CC0;
		l = lbeg+li*dl;
		for(m=0; m<=num-1; ++m)
		{
			cbet = pert(m)+(*mvec[m]).beta;
			cbet = CCI*cbet*(-l);
			am = exp(cbet)*w(m);
			a = a+am;
		}
		p(li) = a.sqabs();
	}
	return p;
}

/* ... write to file */    
void WMM_ModeArray::writeiopower(const WMM_Mode& imode, const WMM_Mode& omode, 
                                Cvector pert, 
                                int numl, double lbeg, double lend,
                                char ext0, char ext1)
{
	Dvector p;
	int i;
	double l, dl;
	char name[13] = "__pow__.xyf";
	FILE *dat;
	if(numl <= 0) wmmarrerror("iopower: numl <= 0");
	switch((*mvec[0]).pol)
	{
		case QTE: name[0] = 't'; name[1] = 'e'; break;
		case QTM: name[0] = 't'; name[1] = 'm'; break;
		case VEC: name[0] = 'v'; name[1] = 'c'; break;
	}
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "P(L = %g -> %g) >> %s\n", lbeg, lend, name);
	p = iopower(imode, omode, pert, numl, lbeg, lend); 
	dat = fopen(name, "w+");
	if(p.nel == 1)
	{
		fprintf(dat, "%g %.10g\n", 0.5*(lbeg+lend), p(0));
		fclose(dat);
		return;
	}
	dl = (lend-lbeg)/(p.nel-1.0);
	for(i=0; i<=p.nel-1; ++i)
	{
		l = lbeg+i*dl;
		fprintf(dat, "%g %.10g\n", l, p(i));
	}
	fclose(dat);
	return;
}

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

/* write single component of the interference field at position z 
   to MATLAB .m file  
   amp: complex amplitudes at z=0
	amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ
   foa: MOD, ORG, SQR, REP, IMP  
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters
   pltype: 'C': contour plot
           'S': surface plot
           'I': intensity image 
           'N': field + mesh only, no plot commands (default) */
void WMM_ModeArray::mfile(Cvector amp, Cvector pert, double z, 
                          Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		          char ext0, char ext1, char pltype) const
{
	FILE *dat;
	char name[13] = "int______.m";
	double rz;
	double max, f, maxsum, minf, maxf; 
	int i;
	rz = floor(z*10.0+0.5)/10.0;
	switch(foa)
	{
		case MOD: name[3] = 'm'; break;
		case ORG: name[3] = 'f'; break;
		case SQR: name[3] = 'q'; break;
		case REP: name[3] = 'r'; break;
		case IMP: name[3] = 'i'; break;
	}
	name[4] = fldchr(cp);
	name[5] = cpchr(cp);
	name[6] = ext0;
	name[7] = ext1;
	name[8] = pltype;
	fprintf(stderr, "%c%c(z=%g) >> %s\n", 
                        fldchr(cp), cpchr(cp), z, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field");
	maxf = 1.0;
	minf = -1.0;
	switch(fldchr(cp))
	{
	case 'E':
		maxsum = 0.0;
		for(i=0; i<=num-1; ++i)
		{
			max = fabs((*mvec[i]).maxE)*amp(i).abs();
			f = fabs((*mvec[i]).minE)*amp(i).abs();
			if(max < f) max = f;
			maxsum += max;
		}
		switch(foa)
		{
		case MOD: case ORG: 
			maxf = maxsum;
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxsum*maxsum;
			minf = 0.0;
			break;
		case REP: case IMP:
			maxf = maxsum;
			minf = -maxsum;
			break;
		}
		break;
	case 'H':
		maxsum = 0.0;
		for(i=0; i<=num-1; ++i)
		{
			max = fabs((*mvec[i]).maxH)*amp(i).abs();
			f = fabs((*mvec[i]).minH)*amp(i).abs();
			if(max < f) max = f;
			maxsum += max;
		}
		switch(foa)
		{
		case MOD: case ORG: 
			maxf = maxsum;
			minf = 0.0;
			break;
		case SQR: 
			maxf = maxsum*maxsum;
			minf = 0.0;
			break;
		case REP: case IMP:
			maxf = maxsum;
			minf = -maxsum;
			break;
		}
		break;
	case 'S':
		maxsum = 0.0;
		for(i=0; i<=num-1; ++i)
		{
			max = fabs((*mvec[i]).maxS)*amp(i).abs()*amp(i).abs();
			f = fabs((*mvec[i]).minS)*amp(i).abs()*amp(i).abs();
			if(max < f) max = f;
			maxsum += sqrt(max);
		}
		switch(foa)
		{
		case MOD: case ORG: case REP: 
			maxf = maxsum*maxsum;
			minf = 0.0;
			break;
		case SQR:  // not reasonable 
			maxf = maxsum*maxsum*maxsum*maxsum;
			minf = 0.0;
			break;
		case IMP:  // not reasonable
			maxf = 1.0;
			minf = 0.0;
		}
		break;
	}
	mlout_geo(dat, (*mvec[0]).wg, minf, maxf);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;
	fld = fieldmat(amp, pert, z, disp, npx, npy, cp, foa);
	mlout_fld(dat, npx, npy, cp, fld);
	name[9] = 0; 
	switch(pltype)
	{
		case 'C': 
			mlout_contour(name, dat, cp, foa); 
			mlout_annotatezpos(dat, rz);
			break;
		case 'S': 
			mlout_surface(name, dat, cp, foa); 
			break;
		case 'I': 
			mlout_image(name, dat, cp, foa); 
			mlout_annotatezpos(dat, rz);
			break;
		default: break;
	}
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}

/* write interference field at position z to MATLAB .m file, fancy style :-)  
   amp: complex amplitudes at z=0
	amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   ext0, ext1: filename id characters */
#define HDIST 1.0e-8
void WMM_ModeArray::fancymfile(Cvector amp, Cvector pert, double z,
                               Rect disp, int npx, int npy, 
                               char ext0, char ext1) const
{
	FILE *dat;
	int np, l, m;
	double x0, x1, xp, y0, y1, yp;
	int numc;
	double max, f, maxsum, minf, maxf; 
	int i;
	Waveguide wg;
	char name[13] = "intfSz__F.m";
	name[6] = ext0;
	name[7] = ext1;
	fprintf(stderr, "Sz(z=%g) >> %s\n", z, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field  :-)");
	name[9] = 0; 
	wg = (*mvec[0]).wg;
	maxsum = 0.0;
	for(i=0; i<=num-1; ++i)
	{
		max = fabs((*mvec[i]).maxS)*amp(i).abs()*amp(i).abs();
		f = fabs((*mvec[i]).minS)*amp(i).abs()*amp(i).abs();
		if(max < f) max = f;
		maxsum += sqrt(max);
	}
	maxf = maxsum*maxsum;
	minf = 0.0;
	mlout_geo(dat, wg, minf, maxf);
	mlout_meshxy(dat, disp, npx, npy);
	Dmatrix fld;
	fld = fieldmat(amp, pert, z, disp, npx, npy, SZ, ORG);
	mlout_fld(dat, npx, npy, SZ, fld);

	Dvector fv;
	numc = 0;
	for(l=0; l<=wg.nx+1; ++l)
	{
		if(l==0) x0 = disp.x0; 	
                else     x0 = wg.hx(l-1);
		if(l==wg.nx+1) x1 = disp.x1; 	
                else           x1 = wg.hx(l);
		for(m=0; m<=wg.ny; ++m)
		{
			if(fabs(wg.n(l,m)-wg.n(l,m+1)) > 1.0e-10)
			{
				yp = wg.hy(m);
				np = (int)(((double) npx)
                                     *(x1-x0)/(disp.x1-disp.x0)); 
				if(np >= 2)
				{
				fv = fieldvec(amp, pert, z, SZ, ORG, np, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp-HDIST, 
                                          x1-HDIST, yp-HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				fv = fieldvec(amp, pert, z, SZ, ORG, np, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST); 
				mlout_sec(dat, 
					  x0+HDIST, yp+HDIST, 
                                          x1-HDIST, yp+HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	for(m=0; m<=wg.ny+1; ++m)
	{
		if(m==0) y0 = disp.y0; 	
                else     y0 = wg.hy(m-1);
		if(m==wg.ny+1) y1 = disp.y1; 	
                else           y1 = wg.hy(m);
		for(l=0; l<=wg.nx; ++l)
		{
			if(fabs(wg.n(l,m)-wg.n(l+1,m)) > 1.0e-10)
			{
				xp = wg.hx(l);
				np = (int)(((double)npy)
                                     *(y1-y0)/(disp.y1-disp.y0)); 
				if(np >= 2)
				{
				fv = fieldvec(amp, pert, z, SZ, ORG, np, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST); 
				mlout_sec(dat, 
					  xp-HDIST, y0+HDIST, 
                                          xp-HDIST, y1-HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				fv = fieldvec(amp, pert, z, SZ, ORG, np, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST);
				mlout_sec(dat, 
					  xp+HDIST, y0+HDIST, 
                                          xp+HDIST, y1-HDIST, 
					  np, SZ, dig10(numc), dig1(numc), fv);
				++numc;
				}
			}
		}
	}
	fv = fieldvec(amp, pert, z, SZ, ORG, 
                       npx, disp.x0, disp.y0, disp.x1, disp.y0); 
	mlout_sec(dat, disp.x0, disp.y0, disp.x1, disp.y0, 
                  npx, SZ, dig10(numc), dig1(numc), fv);
	++numc;
	fv = fieldvec(amp, pert, z, SZ, ORG, 
                       npy, disp.x1, disp.y0, disp.x1, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y0, disp.x1, disp.y1, 
                  npy, SZ, dig10(numc), dig1(numc), fv);
	++numc;
	fv = fieldvec(amp, pert, z, SZ, ORG, 
                       npx, disp.x1, disp.y1, disp.x0, disp.y1); 
	mlout_sec(dat, disp.x1, disp.y1, disp.x0, disp.y1, 
                  npx, SZ, dig10(numc), dig1(numc), fv);
	++numc;
	fv = fieldvec(amp, pert, z, SZ, ORG, 
                       npy, disp.x0, disp.y1, disp.x0, disp.y0); 
	mlout_sec(dat, disp.x0, disp.y1, disp.x0, disp.y0, 
                  npy, SZ, dig10(numc), dig1(numc), fv);
	++numc;

	mlout_fancy(name, dat, SZ, numc);  	
	mlout_print(dat, name, 'p');

	fclose(dat);
	return;
}


/* animate the interference field  
   amp: complex amplitudes at z=0
        amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   cp: EX - HZ, SZ
   foa: MOD, ORG, SQR, REP, IMP 
   disp: output region on the x-y-plane
   npx, npy: number of points in output mesh
   pltype: 'C': contour plot 
           'S': surface plot
           'I': intensity image 
           'F': fancy plot, cp and foa are set to SZ, MOD (default)
   z0, z1, numz: .m files are generated for z-positions 
                  z0+i*(z1-z0)/numz, i=0, ..., numz-1 */
void WMM_ModeArray::movie(Cvector amp, Cvector pert,  
                          Fcomp cp, Afo foa, Rect disp, int npx, int npy, 
		          char pltype,
                          int numz, double z0, double z1) const
{
	int i;
	double z;
	FILE *dat;
	char name[13] = "int___plM.m";
	if(numz >= 100) numz = 100;
	if(pltype == 'F') { cp = SZ; foa = ORG; }
	switch(foa)
	{
		case MOD: name[3] = 'm'; break;
		case ORG: name[3] = 'f'; break;
		case SQR: name[3] = 'q'; break;
		case REP: name[3] = 'r'; break;
		case IMP: name[3] = 'i'; break;
	}
	name[4] = fldchr(cp);
	name[5] = cpchr(cp);
	for(i=0; i<=numz-1; ++i)
	{
		z = z0+i*(z1-z0)/numz; 
		switch(pltype)
		{
		case 'C': case 'S': case 'I':
			mfile(amp, pert, z, cp, foa, 
                              disp, npx, npy, dig10(i), dig1(i), pltype);
			break;
		default:
			fancymfile(amp, pert, z, 
                              disp, npx, npy, dig10(i), dig1(i));
			break;
		}
	}
	fprintf(stderr, ">> %s\n", name);	
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field animation");
	name[8] = pltype;
	name[9] = 0;
	mlout_play(dat, name, numz);
	fclose(dat);
	return;
}


/* interference pattern on the horizontal plane 
   ybeg <= y <= yend, zbeg <= z <= zend at position x
   image plot corresponding to the squareroot of the local intensity (SZ)
   amp: complex amplitudes at z=0
        amplitudes evolve according to amp_j(z) = amp_i(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   npy, npz: number of points in output mesh
   ext0, ext1: filename id characters */
void WMM_ModeArray::prop(Cvector amp, Cvector pert, double x,
                         double ybeg, double yend, int npy,
                         double zbeg, double zend, int npz,
                         char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "propSz__I.m";
	double max, f, maxsum, minf, maxf; 
	double y, dy, z, dz;
	Complex ph;
	int i, j, k, m;
	name[6] = ext0;
	name[7] = ext1;
	fprintf(stderr, "prop(z=%g -> %g) >> %s\n", zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field");
	maxsum = 0.0;
	for(i=0; i<=num-1; ++i)
	{
		max = fabs((*mvec[i]).maxS)*amp(i).abs()*amp(i).abs();
		f = fabs((*mvec[i]).minS)*amp(i).abs()*amp(i).abs();
		if(max < f) max = f;
		maxsum += sqrt(max);
	}
	maxf = maxsum;
	minf = 0.0;
	mlout_geoyz(dat, (*mvec[0]).wg, minf, maxf);
	mlout_meshyz(dat, ybeg, yend, npy, zbeg, zend, npz);

	Dmatrix fld(npy, npz);
	dy = (yend-ybeg)/(npy-1);
	dz = (zend-zbeg)/(npz-1);

	Dmatrix pmprofEx(num, npy);
	Dmatrix pmprofEy(num, npy);
	Dmatrix pmprofHx(num, npy);
	Dmatrix pmprofHy(num, npy);
	for(i=0; i<=npy-1; ++i)
	{
		y = ybeg+i*dy;
		for(m=0; m<=num-1; ++m)
		{
			pmprofEx(m, i) = (*mvec[m]).field(EX, x, y);
			pmprofEy(m, i) = (*mvec[m]).field(EY, x, y);
			pmprofHx(m, i) = (*mvec[m]).field(HX, x, y);
			pmprofHy(m, i) = (*mvec[m]).field(HY, x, y);
		}
	}
	Cvector cbet(num);
	Cvector ma(num);
	for(m=0; m<=num-1; ++m) cbet(m) = pert(m)+(*mvec[m]).beta;
	for(j=0; j<=npz-1; ++j)
	{
		z = zbeg+j*dz;
		for(m=0; m<=num-1; ++m) 
		{
			ph = CCI*cbet(m)*(-z);
			ma(m) = amp(m)*exp(ph);
		}
		for(i=0; i<=npy-1; ++i) 
		{
			y = ybeg+i*dy;
			f = 0.0;
			for(m=0; m<=num-1; ++m) 
			{
				for(k=0; k<=num-1; ++k)
				{
		f += (ma(m).re*ma(k).re+ma(m).im*ma(k).im)
		    *(pmprofEx(m,i)*pmprofHy(k,i)-pmprofEy(m,i)*pmprofHx(k,i));
				}
			}
			fld(i, j) = sqrt(0.5*fabs(f));
		}

	}
	mlout_fld(dat, npy, npz, SZ, fld);
	name[9] = 0; 
	mlout_propimage(name, dat); 
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* Three segment coupler, interference pattern on the horizontal plane 
   ybeg <= y <= yend, zbeg <= z <= zend at elevation x
   image plot corresponding to the squareroot of the local intensity (SZ)
       z < 0: intensity profile of the input modes
   0 < z < l: mode interference pattern
       l < z: intensity profiles of the output modes, 
              scaled by the output amplitudes
(!) this(j) are assumed to be normalized.
   iamp: complex input mode amplitudes
   imodes: input modes, excite the modes this(j) at z=0, 
           belonging to well separated port waveguides 
   iwg: the geometry of the combined input waveguides
   omodes: output modes, belonging to well separated port waveguides 
   owg: the geometry of the combined output waveguides
	iwg and owg are for displaying purposes only
   pert: propagation constant perturbations, cbet(j) = this(j).beta+pert(j)
   npy, npz: number of points in output mesh
   ext0, ext1: filename id characters */
#define ZDIST 1.0e-8
void WMM_ModeArray::ioprop(Cvector iamp, 
                           const WMM_ModeArray& imodes, Waveguide iwg,
                           const WMM_ModeArray& omodes, Waveguide owg,
                           Cvector pert, double x, double l,
                           double ybeg, double yend, int npy,
                           double zbeg, double zend, int npz,
                           char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "propSz__I.m";
	double max, f, maxsum, minf, maxf; 
	int i, j, m, k;
	double y, z, dy, dz;
	Complex a, aj;
	Complex ph;
	double s;
	name[6] = ext0;
	name[7] = ext1;
	if(zbeg >= zend) wmmarrerror("ioprop: invalid z-range");
	fprintf(stderr, "ioprop(z=%g -> %g) >> %s\n", zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "WMM interference field");
	maxsum = 0.0;
	for(i=0; i<=imodes.num-1; ++i)
	{
		max = fabs(imodes(i).maxS)*iamp(i).abs()*iamp(i).abs();
		f = fabs(imodes(i).minS)*iamp(i).abs()*iamp(i).abs();
		if(max < f) max = f;
		maxsum += sqrt(max);
	}
	maxf = maxsum;
	minf = 0.0;
	mlout_iopropgeo(dat, iwg, (*mvec[0]).wg, owg, l, minf, maxf);
	mlout_meshyz(dat, ybeg, yend, npy, zbeg, zend, npz);
	s = 0.0;
	for(i=0; i<=imodes.num-1; ++i)
	{
		fprintf(stderr, "    A^I_%d = %g + i %g,  |A|^2: %g\n", 
                                i, iamp(i).re, iamp(i).im, iamp(i).sqabs());
		s += iamp(i).sqabs();

	}
	fprintf(stderr, "   total input power: %g\n", s); 
	Cvector pamp(num);
	for(i=0; i<=num-1; ++i)
	{
		a = CC0;	
		for(j=0; j<=imodes.num-1; ++j)
		{
			aj = iamp(j)*lscalprod((*mvec[i]), imodes(j));
			a = a+aj;
		}
		pamp(i) = a;
	}
	s=0.0;
	for(i=0; i<=num-1; ++i)
	{
		fprintf(stderr, "    A^II_%d = %g + i %g,  |A|^2: %g\n", 
                                 i, pamp(i).re, pamp(i).im, pamp(i).sqabs());
		s += pamp(i).sqabs();
	}
	fprintf(stderr, "   propagating power: %g\n", s); 
	Cvector oamp(omodes.num);
	for(i=0; i<=omodes.num-1; ++i)
	{
		a = CC0;	
		for(j=0; j<=num-1; ++j)
		{
			aj = pamp(j)*lscalprod(omodes(i), (*mvec[j]));
			ph = pert(j)+(*mvec[j]).beta; 
			ph = CCI*ph*(-l);
			aj = aj*exp(ph);
			a = a+aj;
		}
		oamp(i) = a;
	} 
	s=0.0;
	for(i=0; i<=omodes.num-1; ++i)
	{
		fprintf(stderr, "    A^III_%d = %g + i %g,  |A|^2: %g\n", 
                                 i, oamp(i).re, oamp(i).im, oamp(i).sqabs());
		s += oamp(i).sqabs();
	}
	fprintf(stderr, "   total output power: %g\n", s); 

	Dmatrix fld(npy, npz);
	dy = (yend-ybeg)/(npy-1);
	dz = (zend-zbeg)/(npz-1);

	Dvector improf(npy);
	Dvector omprof(npy);
	for(i=0; i<=npy-1; ++i)
	{
		y = ybeg+i*dy;
		improf(i) = 0.0;
		for(m=0; m<=imodes.num-1; ++m)
		{
			improf(i) += imodes(m).field(SZ, x, y)
		  	              *iamp(m).sqabs();
		}
		omprof(i) = 0.0;
		for(m=0; m<=omodes.num-1; ++m)
		{
			omprof(i) += omodes(m).field(SZ, x, y)
		  	              *oamp(m).sqabs();
		}
	}
	Dmatrix pmprofEx(num, npy);
	Dmatrix pmprofEy(num, npy);
	Dmatrix pmprofHx(num, npy);
	Dmatrix pmprofHy(num, npy);
	for(i=0; i<=npy-1; ++i)
	{
		y = ybeg+i*dy;
		for(m=0; m<=num-1; ++m)
		{
			pmprofEx(m, i) = (*mvec[m]).field(EX, x, y);
			pmprofEy(m, i) = (*mvec[m]).field(EY, x, y);
			pmprofHx(m, i) = (*mvec[m]).field(HX, x, y);
			pmprofHy(m, i) = (*mvec[m]).field(HY, x, y);
		}
	}
	Cvector cbet(num);
	Cvector ma(num);
	for(m=0; m<=num-1; ++m) cbet(m) = pert(m)+(*mvec[m]).beta;
	for(j=0; j<=npz-1; ++j)
	{
		z = zbeg+j*dz;
		if(z<= 0.0)
		{
			for(i=0; i<=npy-1; ++i) fld(i, j) = sqrt(improf(i));
		}
		else
		{
			if(z>=l)
			{
				for(i=0; i<=npy-1; ++i) 
					fld(i, j) = sqrt(omprof(i));
			}
			else
			{
				for(m=0; m<=num-1; ++m) 
				{
					ph = CCI*cbet(m)*(-z);
					ma(m) = pamp(m)*exp(ph);
				}
				for(i=0; i<=npy-1; ++i) 
				{
					y = ybeg+i*dy;
					f = 0.0;
					for(m=0; m<=num-1; ++m) 
					{
						for(k=0; k<=num-1; ++k)
						{
		f += (ma(m).re*ma(k).re+ma(m).im*ma(k).im)
		    *(pmprofEx(m,i)*pmprofHy(k,i)-pmprofEy(m,i)*pmprofHx(k,i));
						}
					}
					fld(i, j) = sqrt(0.5*fabs(f));
				}
			}

		}
	}
	mlout_fld(dat, npy, npz, SZ, fld);
	name[9] = 0; 
	mlout_iopropimage(name, dat); 
	mlout_print(dat, name, 'e');
	fclose(dat);
	return;
}

/* ... as before, but with a single input mode imode only */
void WMM_ModeArray::ioprop(const WMM_Mode& imode, 
                           const WMM_ModeArray& omodes, Waveguide owg,
                           Cvector pert, double x, double l,
                           double ybeg, double yend, int npy,
                           double zbeg, double zend, int npz,
                           char ext0, char ext1) const 
{
	WMM_ModeArray imodes;
	imodes.add(imode);
	Cvector iamp(1);
	iamp(0) = CC1;
	ioprop(iamp, imodes, imode.wg, 
                     omodes, owg, pert, x, l, 
                     ybeg, yend, npy, zbeg, zend, npz, ext0, ext1);
	return;
} 

/* ... as before, with single input and output modes imode and omode */
void WMM_ModeArray::ioprop(const WMM_Mode& imode, 
                           const WMM_Mode& omode, 
                           Cvector pert, double x, double l,
                           double ybeg, double yend, int npy,
                           double zbeg, double zend, int npz,
                           char ext0, char ext1) const
{
	WMM_ModeArray omodes;
	omodes.add(omode);
	ioprop(imode, omodes, omode.wg, pert, x, l, 
               ybeg, yend, npy, zbeg, zend, npz, ext0, ext1);
	return;
}