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

/*
 * slamarr.cpp
 * Arrays of Modes, mode interference evaluation
 */

#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"bepfld.h"
#include"matlvis.h"


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

/* ModeArray: Array of Modes */

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

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

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

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

/* input from FILE dat */
void ModeArray::read(FILE *dat)
{
	if(dat == stderr || dat == stdout) slamarrerror("ModeArray: read");
	if(mvec != NULL) delete[] mvec;
	mvec = NULL;
	num = 0; 
	num = fgetint(dat);
	if(num <= 0) { num = 0; return; }
	mvec = new Mode[num];
	for(int i=0; i<=num-1; ++i)
	{
		mvec[i].read(dat); 
	}
	return;
}

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

/* input from default file */
void ModeArray::read(char ext0, char ext1)
{
	char name[10];
	FILE *dat;
	int i;
	name[0] = 's';
	name[1] = 'l';
	name[2] = 'a';
	name[3] = ext0;
	name[4] = ext1;
	name[5] = '.';
	name[6] = 's';
	name[7] = 'm';
	name[8] = 'a';
	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) %s,  beta: %.10g,  neff: %.10g,  npcB: %.10g\n",
		i, mvec[i].ids, mvec[i].beta, mvec[i].neff, mvec[i].npcB); 
	return;
}

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

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

	return;
}

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

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

/* sort the array by squared propagation constants, highest first */
void ModeArray::sort()
{
	int j, k, maxi;
	double maxb;
	int minn;
	Mode t;
	if(num<=1) return;
	for(j=0; j<=num-2; ++j)
	{
		maxi = j; 
		maxb = mvec[j].betaq;
		minn = mvec[j].nodes();
		for(k=j+1; k<=num-1; ++k)
		{
			if(maxb<mvec[k].betaq-(1.0e-13) ||
			   (fabs(maxb-mvec[k].betaq)<(1.0e-13)
			    && minn > mvec[k].nodes()))
			{
				maxb = mvec[k].betaq;
				minn = mvec[k].nodes();
				maxi = k; 
			}
		}
		t = mvec[j]; mvec[j] = mvec[maxi]; mvec[maxi] = t;
	}
	return;
}

/* compare the present object to ma, return 1 if equal, otherwise return 0 */
int ModeArray::equal(ModeArray& ma)
{
	if(this == &ma) return 1;
	if(num != ma.num) return 0;
	for(int i=0; i<=num-1; ++i) if(mvec[i].equal(ma(i)) != 1) return 0; 
	return 1;
}

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

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

/* bidirectional field superposition at point x, z 
   famp: complex amplitudes of the forward traveling modes at z=0
   bamp: complex amplitudes of the backward traveling modes at z=0
         (z may be negative)
   cp:   EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex ModeArray::field(const Cvector& famp, const Cvector& bamp, 
                         Fcomp cp, double x, double z) const
{
	return field(famp, bamp, cp, x, z, 0.0, 0.0);
}
	
/* bidirectional field superposition at point x, z in [z0, z1],
   famp: complex amplitudes of the forward traveling modes at z=z0
   bamp: complex amplitudes of the backward traveling modes at z=z1
   cp:   EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex ModeArray::field(const Cvector& famp, const Cvector& bamp, 
                         Fcomp cp, double x, double z,
                         double z0, double z1) const
{
	Complex s, ef, eb, f, b;
	Complex ex, ey, ez, hx, hy, hz;
	if(famp.nel < num) slamarrerror("field: array length mismatch");
	if(bamp.nel < num) slamarrerror("field: array length mismatch");
	s = CC0;
	double dzf = z-z0;
	double dzb = z-z1;
	switch(cp)
	{
	case EX: case EY: case EZ: case HX: case HY: case HZ: 
		for(int i=0; i<=num-1; ++i)
		{
			mvec[i].cfield(cp, x, f, b);
			mvec[i].efac(dzf, ef, dzb, eb);
// s += famp(i)*f*ef;  
// s += bamp(i)*b*eb;  
// if(isfinite(ef.re) && isfinite(ef.im)) s += famp(i)*f*ef;  
// if(isfinite(eb.re) && isfinite(eb.im)) s += bamp(i)*b*eb;  
			if(fabs(famp(i).re) > NEGLECT_FLD_CONTRIB ||
			   fabs(famp(i).im) > NEGLECT_FLD_CONTRIB)
				if(std::isfinite(ef.re) && std::isfinite(ef.im)) 
					s += famp(i)*f*ef;  
			if(fabs(bamp(i).re) > NEGLECT_FLD_CONTRIB ||
			   fabs(bamp(i).im) > NEGLECT_FLD_CONTRIB)
				if(std::isfinite(eb.re) && std::isfinite(eb.im)) 
					s += bamp(i)*b*eb;  
		}
		return s;
	case SX:
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		s = ey*hz.conj()-ez*hy.conj();
		return Complex(0.5*s.re, 0.0);
	case SY:
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		s = ez*hx.conj()-ex*hz.conj();
		return Complex(0.5*s.re, 0.0);
	case SZ:
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		s = ex*hy.conj()-ey*hx.conj();
		return Complex(0.5*s.re, 0.0);
	default: slamarrerror("field: mysterious cp");
	}
	return CC0;
}

/* bidirectional field superposition, local field strength at x, z
     id: one of mE, mH, qE, qH, mW
   famp: complex amplitudes of the forward traveling modes at z=0
   bamp: complex amplitudes of the backward traveling modes at z=0
         (z may be negative) */
double ModeArray::lfs(const Cvector& famp, const Cvector& bamp, 
                      FSid id, double x, double z) const
{
	Complex ex, ey, ez, hx, hy, hz;
	double eps_x;
	switch(id)
	{
	case mE: 
		ex = field(famp, bamp, EX, x, z);
		ey = field(famp, bamp, EY, x, z);
		ez = field(famp, bamp, EZ, x, z);
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
	case mH: 
		hx = field(famp, bamp, HX, x, z);
		hy = field(famp, bamp, HY, x, z);
		hz = field(famp, bamp, HZ, x, z);
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
	case qE: 
		ex = field(famp, bamp, EX, x, z);
		ey = field(famp, bamp, EY, x, z);
		ez = field(famp, bamp, EZ, x, z);
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
	case qH: 
		hx = field(famp, bamp, HX, x, z);
		hy = field(famp, bamp, HY, x, z);
		hz = field(famp, bamp, HZ, x, z);
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
	case mW: 
		ex = field(famp, bamp, EX, x, z);
		ey = field(famp, bamp, EY, x, z);
		ez = field(famp, bamp, EZ, x, z);
		hx = field(famp, bamp, HX, x, z);
		hy = field(famp, bamp, HY, x, z);
		hz = field(famp, bamp, HZ, x, z);
		eps_x = mvec[num-1].wg.eps(x);
		return 0.25*(      MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			    +EP0*eps_x*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
		break;
	default: 
		slamarrerror("lfs(), FSid?");
		break;
	}
	return 0.0;
}

/* bidirectional field superposition, local field strength at x, z
   at point x, z in [z0, z1],   
     id: one of mE, mH, qE, qH, mW
   famp: complex amplitudes of the forward traveling modes at z=z0
   bamp: complex amplitudes of the backward traveling modes at z=z1 */
double ModeArray::lfs(const Cvector& famp, const Cvector& bamp, 
                      FSid id, double x, double z, double z0, double z1) const
{
	Complex ex, ey, ez, hx, hy, hz;
	double eps_x;
	switch(id)
	{
	case mE: 
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
	case mH: 
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
	case qE: 
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
	case qH: 
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
	case mW: 
		ex = field(famp, bamp, EX, x, z, z0, z1);
		ey = field(famp, bamp, EY, x, z, z0, z1);
		ez = field(famp, bamp, EZ, x, z, z0, z1);
		hx = field(famp, bamp, HX, x, z, z0, z1);
		hy = field(famp, bamp, HY, x, z, z0, z1);
		hz = field(famp, bamp, HZ, x, z, z0, z1);
		eps_x = mvec[num-1].wg.eps(x);
		return 0.25*(      MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			    +EP0*eps_x*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
		break;
	default: 
		slamarrerror("lfs(z0, z1), FSid?");
		break;
	}
	return 0.0;
}

/* bidirectional field superposition, local energy density at point x, z
   famp: complex amplitudes of the forward traveling modes at z=0
   bamp: complex amplitudes of the backward traveling modes at z=0
         (z may be negative) */
double ModeArray::endens(const Cvector& famp, const Cvector& bamp, 
                         double x, double z) const
{
	return lfs(famp, bamp, mW, x, z);
}

/* bidirectional field superposition, 
   local energy density at point x, z in [z0, z1],
   famp: complex amplitudes of the forward traveling modes at z=z0
   bamp: complex amplitudes of the backward traveling modes at z=z1 */
double ModeArray::endens(const Cvector& famp, const Cvector& bamp, 
                         double x, double z, double z0, double z1) const
{
	return lfs(famp, bamp, mW, x, z, z0, z1);
}

/* the six electric and magnetic components at x, z */
void ModeArray::emfield(const Cvector& famp, const Cvector& bamp,
                        double x, double z, double z0, double z1,
                        Complex& ex, Complex& ey, Complex& ez,
                        Complex& hx, Complex& hy, Complex& hz) const
{
	Complex s, ef, eb, f, b;
	if(famp.nel < num) slamarrerror("field: array length mismatch");
	if(bamp.nel < num) slamarrerror("field: array length mismatch");

	ex = ey = ez = hx = hy = hz = CC0;
	double dzf = z-z0;
	double dzb = z-z1;
	Complex af, ab;
	bool icf, icb;
	for(int i=0; i<=num-1; ++i)
	{
		af = famp(i);
		ab = bamp(i);
		mvec[i].efac(dzf, ef, dzb, eb);

		icf = true;
		if(fabs(af.re) < NEGLECT_FLD_CONTRIB && 
		   fabs(af.im) < NEGLECT_FLD_CONTRIB) icf = false;
		if(!std::isfinite(ef.re)) icf = false;  
		if(!std::isfinite(ef.im)) icf = false;  
		icb = true;
		if(fabs(ab.re) < NEGLECT_FLD_CONTRIB && 
		   fabs(ab.im) < NEGLECT_FLD_CONTRIB) icb = false;
		if(!std::isfinite(eb.re)) icb = false;  
		if(!std::isfinite(eb.im)) icb = false;  

		if(icf==false && icb==false) return;

		af *= ef;
		ab *= eb;

		mvec[i].cfield(EX, x, f, b);
		if(icf) ex += af*f;  
		if(icb) ex += ab*b;  

		mvec[i].cfield(EY, x, f, b);
		if(icf) ey += af*f;  
		if(icb) ey += ab*b;  

		mvec[i].cfield(EZ, x, f, b);
		if(icf) ez += af*f;  
		if(icb) ez += ab*b;  

		mvec[i].cfield(HX, x, f, b);
		if(icf) hx += af*f;  
		if(icb) hx += ab*b;  

		mvec[i].cfield(HY, x, f, b);
		if(icf) hy += af*f;  
		if(icb) hy += ab*b;  

		mvec[i].cfield(HZ, x, f, b);
		if(icf) hz += af*f;  
		if(icb) hz += ab*b;  
	}
	return;
}

void ModeArray::emfield(const Cvector& famp, const Cvector& bamp,
                        double x, double z,
                        Complex& ex, Complex& ey, Complex& ez,
                        Complex& hx, Complex& hy, Complex& hz) const
{
	emfield(famp, bamp, x, z, 0.0, 0.0, ex, ey, ez, hx, hy, hz);
	return;
}

/* field superposition at point (x, z), 
   amp:  complex amplitudes of forward traveling waves at z=0,
	 amplitudes evolve according to amp_j(z) = amp_j(0)*exp(-i cbet_j*z)
   pert: propagation constant perturbations, with
         cbet_j = this(j).cbeta(FORW)+pert(j)
   cp:   EX, EY, EZ, HX, HY, HZ, SX, SY, SZ */
Complex ModeArray::pfield(const Cvector& amp, const Cvector& pert, Fcomp cp, 
                          double x, double z) const

{
	Complex s, e, f;
	Complex ex, ey, ez, hx, hy, hz;
	if( amp.nel < num) slamarrerror("pfield: array length mismatch");
	if(pert.nel < num) slamarrerror("pfield: array length mismatch");
	s = CC0;
	switch(cp)
	{
	case EX: case EY: case EZ: case HX: case HY: case HZ: 
		for(int i=0; i<=num-1; ++i)
		{
			e = mvec[i].efac(FORW, pert(i), z);
			f = mvec[i].cfield(cp, FORW, x);
			s += amp(i)*f*e;  
		}
		return s;
	case SX:
		ey = pfield(amp, pert, EY, x, z);
		ez = pfield(amp, pert, EZ, x, z);
		hy = pfield(amp, pert, HY, x, z);
		hz = pfield(amp, pert, HZ, x, z);
		s = ey*hz.conj()-ez*hy.conj();
		return Complex(0.5*s.re, 0.0);
	case SY:
		ex = pfield(amp, pert, EX, x, z);
		ez = pfield(amp, pert, EZ, x, z);
		hx = pfield(amp, pert, HX, x, z);
		hz = pfield(amp, pert, HZ, x, z);
		s = ez*hx.conj()-ex*hz.conj();
		return Complex(0.5*s.re, 0.0);
	case SZ:
		ex = pfield(amp, pert, EX, x, z);
		ey = pfield(amp, pert, EY, x, z);
		hx = pfield(amp, pert, HX, x, z);
		hy = pfield(amp, pert, HY, x, z);
		s = ex*hy.conj()-ey*hx.conj();
		return Complex(0.5*s.re, 0.0);
	default: slamarrerror("pfield: mysterious cp");
	}
	return CC0;
}

/* assuming that all modes belong to the same mode spectrum:
   check orthogonality & normalization */
double ModeArray::orthonormality()
{
	double devn = 0.0;
	double devo = 0.0;
	double o; 
	Mode mj;
	for(int j=0; j<=num-1; ++j)
	{
		mj = mvec[j];
		if(mj.ppt == PROPAG)
		{
			o = (mj.overlap(FORW, mj, FORW)).re;
			if(fabs(1-o) > devn) devn = fabs(1-o);
		}
		else 
		{
			o = (mj.overlap(FORW, mj, BACK)).abs();
			if(fabs(1-o) > devn) devn = fabs(1-o);
		}
	}
	for(int j=0; j<=num-1; ++j)
	{
		mj = mvec[j];
		for(int k=0; k<=num-1; ++k)
		{
			o = (mj.overlap(FORW, mvec[k], FORW)).abs();
			if(j != k)
			{
				if(o > devo) devo = o;
			}
		}
	}
	fprintf(stderr, "DevNorm: %g, DevOrthog: %g\n", devn, devo);
	if(devo > devn) return devo;
	return devn;
}

/* related to the classification of FB modes:
   remove all non propagating modes from the array */
void ModeArray::restrict2prop()
{	
	int done = 0;
	while(done == 0 && num >= 1)
	{
		done=1;
		int r=-1;
		for(int j=0; j<=num-1; ++j)
		{
			if((mvec[j]).ppt==EVANESC)
			{
				done = 0;
				r = j; 
			}
		}
		if(done == 0) remove(r);
	}
	return;
}


/* - 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 ModeArray::pweight(const Mode& imode, int m,
                          const Mode& omode) 
{
	Interval i(-AWAY, AWAY);
	double om, mi;
	om = 0.25*(omode.integrate(EX, (mvec[m]), HY, i)
	          -omode.integrate(EY, (mvec[m]), HX, i)
	          +(mvec[m]).integrate(EX, omode, HY, i)
	          -(mvec[m]).integrate(EY, omode, HX, i));
	mi = 0.25*((mvec[m]).integrate(EX, imode, HY, i)
	          -(mvec[m]).integrate(EY, imode, HX, i)
	          +imode.integrate(EX, (mvec[m]), HY, i)
	          -imode.integrate(EY, (mvec[m]), HX, i));
	return om*mi;
}
   
/* single relative power level for device length l */    
double ModeArray::iopower(const Mode& imode, const Mode& omode, 
                          const 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 ModeArray::iopower(const Mode& imode, const Mode& omode, 
                           const 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) slamarrerror("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 ModeArray::writeiopower(const Mode& imode, const Mode& omode, 
                             const Cvector& pert, 
                             int numl, double lbeg, double lend,
                             char ext0, char ext1)
{
	Dvector p;
	int i;
	double l, dl;
	char name[13] = "t_pow__.xyf";
	FILE *dat;
	if(numl <= 0) slamarrerror("iopower: numl <= 0");
	name[1] = polchr(mvec[0].pol);
	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 ---------------------------------------- */

/* mode interference plots, image of component cp
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
   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 ModeArray::plot(const Cvector& amp_f, const Cvector& amp_b, 
                     Fcomp cp, Afo foa,
                     double xbeg, double xend, double zbeg, double zend, 
                     int npx, int npz, char ext0, char ext1) const
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
	fld.plot(cp, foa, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

/* image of the relative phase
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   cp: one of EX, EY, EZ, HX, HY, HZ
   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 ModeArray::phasemap(const Cvector& amp_f, const Cvector& amp_b, 
                         Fcomp cp,
                         double xbeg, double xend, double zbeg, double zend,
                         int npx, int npz, char ext0, char ext1) const
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
	fld.phasemap(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

/* electromagnetic energy density image
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   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 ModeArray::plot(const Cvector& amp_f, const Cvector& amp_b, 
                     double xbeg, double xend, double zbeg, double zend,
                     int npx, int npz, char ext0, char ext1) const
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
	fld.plot(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

/* mode interference plots, fancy style surface ;-)
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
   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 ModeArray::fplot(const Cvector& amp_f, const Cvector& amp_b, 
                      Fcomp cp,
                      double xbeg, double xend, double zbeg, double zend, 
                      int npx, int npz, char ext0, char ext1) const
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
	fld.fplot(cp, xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
	return;
}

/* animation of the light propagation along the waveguide,
   the frames show images of EY (TE) resp. HY (TM), at
   ntfr equidistant times over one time period
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ntfr: number of frames per time period
   ext0, ext1: filename id characters for the m.file */
void ModeArray::movie(const Cvector& amp_f, const Cvector& amp_b,
                      double xbeg, double xend, double zbeg, double zend,
                      int npx, int npz, int ntfr, char ext0, char ext1)
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
	fld.movie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1);
	return;
}

/* animation of the light propagation along the waveguide, fancy style ;-)
   the frames show surfaces of EY (TE) resp. HY (TM), at
   ntfr equidistant times over one time period
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window   
   npx, npz: number of plot points in the x and z directions
   ntfr: number of frames per time period
   ext0, ext1: filename id characters for the m.file */
void ModeArray::fmovie(const Cvector& amp_f, const Cvector& amp_b,
                       double xbeg, double xend, double zbeg, double zend,
                       int npx, int npz, int ntfr, char ext0, char ext1)
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
	fld.fmovie(xbeg, xend, zbeg, zend, npx, npz, ntfr, ext0, ext1);
	return;
}

/* export full field data ("all") into a viewer m-file
   amp_f: amplitudes of the forward propagating modes at z=0 
   amp_b: amplitudes of the backward propagating modes at z=0
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters */
void ModeArray::viewer(const Cvector& amp_f, const Cvector& amp_b,
                       double xbeg, double xend, double zbeg, double zend,
                       int npx, int npz, char ext0, char ext1) const
{
	SegWgStruct st(-1);
	st(0) = mvec[num-1].wg;
	BepField fld(st);
	fld(0) = *this;
	fld.setf(0, SBRIGHT, amp_f);
	fld.setf(0, SBLEFT, amp_f);
	fld.setb(0, SBRIGHT, amp_b);
	fld.setb(0, SBLEFT, amp_b);
	fld.check();
        fld.viewer(xbeg, xend, zbeg, zend, npx, npz, ext0, ext1);
        return;
}

/* - for the QUEP implementation ---------------------------------------- */

/* overlaps of perpendiculary propagating modes */
Cmatrix ModeArray::crossoverlap(Propdir hdir, 
                                const ModeArray& mv, Propdir vdir,
                                double zpos, Interval i) const
{
fprintf(stderr, "marr: cro: used.\n");
	Cmatrix g(num, mv.num);
	for(int r=0; r<=num-1; ++r)
	{
		for(int m=0; m<=mv.num-1; ++m)
		g(r, m) = mvec[r].crossoverlap(hdir, mv(m), vdir, zpos, i);
	}
	return g;
}

void ModeArray::crossoverlap(const ModeArray& mv,
                             double zpos, Interval i,
                             Cmatrix& cff, Cmatrix& cfb, 
                             Cmatrix& cbf, Cmatrix& cbb) const
{
	cff = Cmatrix(num, mv.num);
	cfb = Cmatrix(num, mv.num);
	cbf = Cmatrix(num, mv.num);
	cbb = Cmatrix(num, mv.num);
	for(int r=0; r<=num-1; ++r)
	{
		for(int m=0; m<=mv.num-1; ++m)
		{
			mvec[r].crossoverlap(mv(m), zpos, i,
		                             cff(r, m), cfb(r, m), 
		                             cbf(r, m), cbb(r, m));
		}
	}
	return;
}

/* - for the vQUEP implementation --------------------------------------- */

/* prepare vTE / vTM mode representation, for given ky, for all modes */
void ModeArray::vectorify(Complex ky)
{
	for(int r=0; r<=num-1; ++r) mvec[r].vectorify(ky);
	return;	
}