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

/*
 * eim.cpp
 * Procedures related to an effective-index-like dimensionality reduction
 * 2D -> 1D for integrated optical scattering problems
 */

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

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

/* ------------------------------------------------------------------------
   field container, including the solver procedure 
   ------------------------------------------------------------------------ */

/* initialize */
EimField::EimField(SegWgStruct s, Polarization p, double wl)
{
	str = s;
	pol = p;
	lambda = wl;
	k0 = val_k0(lambda);
	invommu0 = val_invommu0(lambda);
	invomep0 = val_invomep0(lambda);
}

/* the vEIM solver */
void EimField::solve(int variant)
{
	fprintf(stderr, "\n");
	fprintf(stderr, "------------- vEIM ------------------------------------------------- '09 ---\n");
fprintf(stderr, "T%c  lambda: %.10g  k0: %g", polCHR(pol), lambda, k0);
if(str.special()) fprintf(stderr, "  (!)\n");
else           fprintf(stderr, "\n");
	fprintf(stderr, "----------------------------------------------------------------------------\n");
	Circuit strc = str.circuit();
fprintf(stderr, "    Nx: %d\n", strc.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=strc.nx; ++j) fprintf(stderr, "%6.10g  ", strc.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    Nz: %d\n", strc.nz);
fprintf(stderr, "    hz: ");
for(int j=0; j<=strc.nz; ++j) fprintf(stderr, "%6.10g  ", strc.hz(j));
fprintf(stderr, "\n");
if(strc.special) fprintf(stderr, "   eps:\n");
else           fprintf(stderr, "     n:\n");
for(int i=strc.nx+1; i>=0; --i) 
{
fprintf(stderr,"layer %d:  ", i);
for(int j=0; j<=strc.nz+1; ++j) fprintf(stderr, "%6.10g  ", strc.n(i, j));
fprintf(stderr,"\n");
}
	fprintf(stderr, "--- reference --------------------------------------------------------------\n");

	Waveguide rs = str(0);
	ModeArray ma;
	modeanalysis(rs, pol, ma);
	if(ma.num < 1) eimerror("no guided mode in the reference slice"); 
	vp = ma(0);
	if(variant == 0)
		fprintf(stderr, "--- effective permittivity -------------------------------------------------\n");
	else
		fprintf(stderr, "--- effective permittivity - (!) -------------------------------------------\n");
	epseff = Waveguide(str.nz);
	epseff.hx = str.hz;
	epseff.lambda = str(0).lambda;
	epseff.special = 1;
	epseff_b = epseff;

	if(pol == TE)
	{
		for(int s=0; s<=strc.nz+1; ++s)
		{
			double na = 0.0;
			double no = 0.0;
			for(int l=0; l<=strc.nx+1; ++l)
			{
				double x0, x1;
				if(l==0) x0 = -AWAY; 
				else     x0 = strc.hx(l-1);
				if(l==strc.nx+1) x1 = AWAY; 
				else             x1 = strc.hx(l);
				double i = vp.integrate(EY, EY, 
				                        Interval(x0, x1));
				no += i;
				na += i*(strc.n(l, s)*strc.n(l, s) 
					 - rs.eps(0.5*(x1+x0)));
			}
			epseff.n(s) = vp.neff*vp.neff+na/no;
		}
		epseff_b.n.init(1.0);
	}
	else
	{
if(variant == 0)
{
		double na = 0.0;
		for(int l=0; l<=rs.nx+1; ++l)
		{
			double x0, x1;
			if(l==0) x0 = -AWAY; 
			else     x0 = rs.hx(l-1);
			if(l==rs.nx+1) x1 = AWAY; 
			else           x1 = rs.hx(l);
			double i = vp.integrate(HY, HY, 
						Interval(x0, x1));
			na += i/(rs.n(l)*rs.n(l)); 
		}
		for(int s=0; s<=strc.nz+1; ++s)
		{
			double nb = 0.0;
			for(int l=0; l<=strc.nx+1; ++l)
			{
				double x0, x1;
				if(l==0) x0 = -AWAY; 
				else     x0 = strc.hx(l-1);
				if(l==strc.nx+1) x1 = AWAY; 
				else             x1 = strc.hx(l);
				double i = vp.integrate(HY, HY, 
				                        Interval(x0, x1));
				nb += i/(strc.n(l, s)*strc.n(l, s)); 
			}
			epseff_b.n(s) = na/nb;

			double nd = 0.0;
			for(int l=0; l<=strc.nx+1; ++l)
			{
				double x0, x1;
				if(l==0) x0 = -AWAY; 
				else     x0 = strc.hx(l-1);
				if(l==strc.nx+1) x1 = AWAY; 
				else             x1 = strc.hx(l);
				double i = vp.integrate(EZ, EZ, 
				                        Interval(x0, x1));
				double er = rs.eps(0.5*(x0+x1));
				double e  = strc.n(l, s)*strc.n(l, s); 
				nd += i*(1.0/er - 1.0/e)*er*er/
				      (invomep0*invomep0); 
			}
			epseff.n(s) = vp.neff*vp.neff*na/nb+nd/nb/k0/k0;
		}
}
else
{
		for(int s=0; s<=strc.nz+1; ++s)
		{
			double na = 0.0;
			double nb = 0.0;
			double nc = 0.0;
			double nd = 0.0;
			double ne = 0.0;
			double nf = 0.0;
			double ng = 0.0;
			double nh = 0.0;
			for(int l=0; l<=strc.nx+1; ++l)
			{
				double x0, x1;
				if(l==0) x0 = -AWAY; 
				else     x0 = strc.hx(l-1);
				if(l==strc.nx+1) x1 = AWAY; 
				else             x1 = strc.hx(l);
				double er = rs.eps(0.5*(x0+x1));
				double e  = strc.n(l, s)*strc.n(l, s); 
				na += er * 
				      vp.integrate(EZ, EZ, Interval(x0, x1));
				nb += e  * 
				      vp.integrate(EZ, EZ, Interval(x0, x1));
				nc +=      
				      vp.integrate(HY, HY, Interval(x0, x1));
				nd += 1.0/er * 
				      vp.integrate(HY, HY, Interval(x0, x1));
				ne += (e-er) * 
				      vp.integrate(EZ, EZ, Interval(x0, x1));
				nf += e * 
				      vp.integrate(EZ, EZ, Interval(x0, x1));
				ng += e * 
				      vp.integrate(EX, EX, Interval(x0, x1));
				nh += er * 
				      vp.integrate(EX, EX, Interval(x0, x1));
			}
			epseff.n(s) = (vp.neff*vp.neff*na/nb+nc*ne/nd/nf)*ng/nh;
			epseff_b.n(s) = ng/nh;
		}
		
}
	}

	for(int s=0; s<=strc.nz+1; ++s)
	{
		fprintf(stderr, "(%d)  %g  %g\n", 
			s, epseff.n(s), epseff_b.n(s));
/* check epseff 
double n1 = 0, d = 0, n2 = 0;	
double x0 = strc.hx(0)-7.0;
double x1 = strc.hx(strc.nx)+7.0;
double dx = 0.0001;
double er, e, del, chi, dchi;
for(double x=x0; x<=x1; x+=dx) 
{
	chi = vp.field(HY, x);
	dchi = (vp.field(HY, x+dx/20.0)
	       -vp.field(HY, x-dx/20.0))/(dx/10.0);
	er = rs.eps(x);
	e  = str(s).eps(x);
	del = 1.0/(1.0/er-1.0/e);
	n1 += chi*chi/er;
	d  += chi*chi/e;
	n2 += dchi*dchi/del;
}
double r = (vp.beta/k0)*(vp.beta/k0)*n1/d+1.0/(k0*k0)*n2/d;
fprintf(stderr, "(%d)  %g  %g\n", s, r, n1/d);
*/

	}

	fprintf(stderr, "----------------------------------------------------------------------------\n");

	// ---- transfer matrix solver, 1D multilayer ----

	double in = 1.0;
	int num;
	Complex df;
	Cvector dum(2);
	double err;

	num = epseff.nx;

	Cvector u(num+2);
	Cvector r(num+2);
	Dvector xu(num+2);
	Dvector xr(num+2);
	Dvector gamma(num+2);
	Cvector cg(num+2);
	Ivector s(num+2);
	Dvector eta(num+2);
	
	Cvector fca(num+2);
	Cvector fcb(num+2);
	Cvector fcc(num+2);
	Cvector fcd(num+2);
	Cvector fcA(num+2);
	Cvector fcB(num+2);
	Cvector fcC(num+2);
	Cvector fcD(num+2);
	Cvector fcR(num+2);

	for(int j=0; j<=num+1; ++j)
	{
		double g = k0*k0*epseff.eps(j);
		gamma(j)=sqrt(fabs(g));
		if(g >= 0.0)
		{
			s(j) = 1;
			cg(j) = Complex(gamma(j), 0.0);
			if(j==0) xu(j) = xr(j) = epseff.hx(0);
			else
			{
				if(j==num+1) xu(j)=xr(j)=epseff.hx(num);
				else xr(j)=xu(j)=0.5*(epseff.hx(j-1)+epseff.hx(j)); 
			}
		}
		else
		{
			s(j) = -1;
			cg(j) = Complex(0.0, -gamma(j));
			if(j >= 1  ) xu(j)=epseff.hx(j-1); else xu(j)=epseff.hx(0);
			if(j <= num) xr(j)=epseff.hx(j);   else xr(j)=epseff.hx(num);
		}
		if(pol == TE) eta(j)=1.0;
		else          eta(j)=1.0/epseff_b.eps(j);
	}	
	for(int j=0; j<=num; ++j)
	{
		double dxu, dxr;
		Complex egu, egr;
		dxu = epseff.hx(j)-xu(j);
		dxr = epseff.hx(j)-xr(j);
		egu = exp(CCI*cg(j)*(-dxu));
		egr = exp(CCI*cg(j)*( dxr));
		fca(j) = egu; 
		fcb(j) = egr; 
		fcc(j) = CCI*cg(j)*egu*(-eta(j));
		fcd(j) = CCI*cg(j)*egr*( eta(j));

		dxu = epseff.hx(j)-xu(j+1);
		dxr = epseff.hx(j)-xr(j+1);
		egu = exp(CCI*cg(j+1)*(-dxu));
		egr = exp(CCI*cg(j+1)*( dxr));
		fcA(j) = egu; 
		fcB(j) = egr; 
		fcC(j) = CCI*cg(j+1)*egu*(-eta(j+1));
		fcD(j) = CCI*cg(j+1)*egr*( eta(j+1));
	}	

	fcR(num) = (fcc(num)/fcC(num)-fca(num)/fcA(num))/
	           (fcb(num)/fcA(num)-fcd(num)/fcC(num));
	for(int j=num-1; j>=0; --j)
	{	
		fcR(j) = ( fcc(j)/(fcC(j)+fcD(j)*fcR(j+1))
                          -fca(j)/(fcA(j)+fcB(j)*fcR(j+1)))/
		         ( fcb(j)/(fcA(j)+fcB(j)*fcR(j+1))
                          -fcd(j)/(fcC(j)+fcD(j)*fcR(j+1)));
	}

	if(pol == TE) u(0) = Complex(sqrt(2.0*fabs(in)
	                                  /gamma(0)/invommu0));
	else          u(0) = Complex(sqrt(2.0*fabs(in)*epseff.eps(0)
	                                  /gamma(0)/invomep0));
	r(0) = fcR(0)*u(0);
	for(int j=1; j<=num; ++j)
	{
		u(j) = (fca(j-1)*u(j-1)+fcb(j-1)*r(j-1))/
		       (fcA(j-1)+fcB(j-1)*fcR(j));
		r(j) = fcR(j)*u(j);
	}
	u(num+1) = (fca(num)*u(num)+fcb(num)*r(num))/fcA(num);
	r(num+1) = CC0;

	err = 0.0;
	for(int j=0; j<=num; ++j)
	{
	df = 
	(u(j)*exp(CCI*cg(j)*(-(epseff.hx(j)-xu(j))))
	+r(j)*exp(CCI*cg(j)*( (epseff.hx(j)-xr(j)))))-
	(u(j+1)*exp(CCI*cg(j+1)*(-(epseff.hx(j)-xu(j+1))))
	+r(j+1)*exp(CCI*cg(j+1)*( (epseff.hx(j)-xr(j+1)))));
	err += df.abs();
	df = 
	(u(j)*exp(CCI*cg(j)*(-(epseff.hx(j)-xu(j))))*CCI*cg(j)*(-eta(j))
	+r(j)*exp(CCI*cg(j)*( (epseff.hx(j)-xr(j))))*CCI*cg(j)*( eta(j)))-
	(u(j+1)*exp(CCI*cg(j+1)*(-(epseff.hx(j)-xu(j+1))))*CCI*cg(j+1)*(-eta(j+1))
	+r(j+1)*exp(CCI*cg(j+1)*( (epseff.hx(j)-xr(j+1))))*CCI*cg(j+1)*( eta(j+1)));
	err += df.abs();
	}
	if(err/(u(0).abs()) > 1.0e-8) eimerror("accuracy");

	Mode mfr, mfi;
	mfr.wg = epseff;
	mfi.wg = epseff;
	mfr.pol = pol;
	mfi.pol = pol;
	mfr.k0 = k0;
	mfi.k0 = k0;
	mfr.invommu0 = invommu0;
	mfi.invommu0 = invommu0;
	mfr.invomep0 = invomep0;
	mfi.invomep0 = invomep0;
	mfr.beta = 0.0;
	mfi.beta = 0.0;
	mfr.betaq = 0.0;
	mfi.betaq = 0.0;
	mfr.neff = 0.0;
	mfi.neff = 0.0;
	mfr.npcB = 0.0;
	mfi.npcB = 0.0;
	mfr.ppt = PROPAG;
	mfi.ppt = PROPAG;
	mfr.bdt_b = OPEN;
	mfi.bdt_b = OPEN;
	mfr.bp_b = -AWAY;
	mfi.bp_b = -AWAY;
	mfr.bdt_t = OPEN;
	mfi.bdt_t = OPEN;
	mfr.bp_t = AWAY;
	mfi.bp_t = AWAY;

	mfr.N = num+2;
	mfi.N = num+2;
	mfr.piece = new SmPiece[num+2];
	mfi.piece = new SmPiece[num+2];
	for(int l=0; l<=num+1; ++l)
	{
		mfr.piece[l].pol = pol;
		mfi.piece[l].pol = pol;
		Interval i = epseff.layer(l);
		mfr.piece[l].bif = i.x0;
		mfi.piece[l].bif = i.x0;
		mfr.piece[l].tif = i.x1;
		mfi.piece[l].tif = i.x1;
		mfr.piece[l].k0 = k0;
		mfi.piece[l].k0 = k0;
		mfr.piece[l].k0nq = k0*k0*epseff.eps(l);
		mfi.piece[l].k0nq = k0*k0*epseff.eps(l);
		mfr.piece[l].inve = 1.0/epseff.eps(l);
		mfi.piece[l].inve = 1.0/epseff.eps(l);
		
		mfr.piece[l].g = gamma(l);
		mfi.piece[l].g = gamma(l);
		mfr.piece[l].xa = xr(l);
		mfi.piece[l].xa = xr(l);
		mfr.piece[l].xb = xu(l);
		mfi.piece[l].xb = xu(l);

		double gq = -k0*k0*epseff.eps(l);
		if(gq > 0.0) 
		{
			mfr.piece[l].bhv = HYP;
			mfi.piece[l].bhv = HYP;
		}
		else	
		{
			mfr.piece[l].bhv = OSC;
			mfi.piece[l].bhv = OSC;
		}
	}
	for(int j=0; j<=num+1; ++j)  
	{
		if(s(j) == 1)
		{
			mfr.piece[j].A = -(r(j)-u(j)).im;
			mfr.piece[j].B =  (r(j)+u(j)).re;
			mfi.piece[j].A =  (r(j)-u(j)).re;
			mfi.piece[j].B =  (r(j)+u(j)).im;
		}
		else
		{
			mfr.piece[j].A = r(j).re;
			mfr.piece[j].B = u(j).re;
			mfi.piece[j].A = r(j).im;
			mfi.piece[j].B = u(j).im;
		}
	}
	mfr.special = -1;
	mfi.special = -1;
	mfr.ids[0] = 'p'; mfr.ids[1] = 'w';
	mfr.ids[2] = 'R'; mfr.ids[3] = 'e';
	mfr.ids[4] = 0; mfr.ids[5] = 0; mfr.ids[6] = 0; mfr.ids[7] = 0;
	mfi.ids[0] = 'p'; mfi.ids[1] = 'w';
	mfi.ids[2] = 'I'; mfi.ids[3] = 'm';
	mfi.ids[4] = 0; mfi.ids[5] = 0; mfi.ids[6] = 0; mfi.ids[7] = 0;
	mfr.setfieldmap();
	mfi.setfieldmap();
	hp.clear();
	hp.add(mfr);
	hp.add(mfi);
	amp = Cvector(2);
	amp(0) = CC1;
	amp(1) = CCI;
	dum.init(CC0);

	if(pol == TE) ref = 0.5*gamma(0)*invommu0*r(0).sqabs();
	else ref = 0.5*gamma(0)*invomep0*r(0).sqabs()/epseff.eps(0);

	if(pol == TE)
	{
		if(s(num+1) > 0) trans = 0.5*gamma(num+1)
		                 *invommu0*u(num+1).sqabs();
		else trans = 0.0;
	}
	else
	{
		if(s(num+1) > 0) trans = 0.5*gamma(num+1)/epseff.eps(num+1)
		                 *invomep0*u(num+1).sqabs();
		else trans = 0.0;
	}
	if(fabs(in-ref-trans) > 1.0e-10) eimerror("energy conservation");

	fprintf(stderr, "T  = %g\n", trans);
	fprintf(stderr, "R  = %g\n", ref);
	fprintf(stderr, "Lo = %g\n", 1.0-ref-trans);
	fprintf(stderr, "----------------------------------------------------------------------------\n");

/* check wave eqn 
double z0 = str.hz(0)-5.0;
double z1 = str.hz(str.nz)+5.0;
double dz = 0.1;
double ddz = 0.0001;
double re; 
Complex psi, psim, psip, ddpsi, errr;
for(double z=z0; z<=z1; z+=dz) 
{
	psi  = amp(0)*hp(0).field(HY,z)+amp(1)*hp(1).field(HY,z);
	psim = amp(0)*hp(0).field(HY,z-ddz)+amp(1)*hp(1).field(HY,z-ddz);
	psip = amp(0)*hp(0).field(HY,z+ddz)+amp(1)*hp(1).field(HY,z+ddz);
	ddpsi = (psim+psip-2.0*psi)/ddz/ddz;
	re = epseff.eps(z);
	errr = ddpsi+k0*k0*re*psi;
	fprintf(stderr, "z=%g  err=%g\n", z, errr.abs());
}
for(int s=0; s<=str.nz; ++s) 
{ 
	double z = str.hz(s);
	psi  = amp(0)*hp(0).field(HY,z)+amp(1)*hp(1).field(HY,z);
	psim = amp(0)*hp(0).field(HY,z-ddz)+amp(1)*hp(1).field(HY,z-ddz);
	psip = amp(0)*hp(0).field(HY,z+ddz)+amp(1)*hp(1).field(HY,z+ddz);
	Complex dpsim, dpsip;
	double eem = epseff_b.eps(z-ddz);
	double eep = epseff_b.eps(z+ddz);
	dpsim = (psi-psim)/ddz;
	dpsip = (psip-psi)/ddz;
	errr = dpsip/eep-dpsim/eem; 
	fprintf(stderr, "if %d  err=%g\n", s, errr.abs());
}
exit(1);
*/

	return;
}


/* the resulting field */
Complex EimField::field(Fcomp cp, double x, double z) const
{
	Complex h, v, c;
	if(pol == TE)
	{	
		switch(cp)
		{
		case EX: return CC0;
			break;
		case EY: 
			h = amp(0)*hp(0).field(EY,z)+amp(1)*hp(1).field(EY,z);
			return h*vp.field(EY, x);
			break;
		case EZ: return CC0;
			break;
		case HX: 
			h = (amp(0)*hp(0).field(HZ,z)+amp(1)*hp(1).field(HZ,z))
			    /invommu0;
			v = vp.field(EY, x);
			c = CCmI*invommu0;
			return c*h*v;
			break;
		case HY: return CC0;
			break;
		case HZ:
			h = amp(0)*hp(0).field(EY,z)+amp(1)*hp(1).field(EY,z);
			v = vp.field(HZ, x)/invommu0;
			c = CCI*invommu0;
			return c*h*v;
			break;
		default: eimerror("field: cp not implemented"); 
			break;
		}
		return CC0;
	}
	switch(cp)
	{
	case EX: 
		h = (amp(0)*hp(0).field(EZ,z)+amp(1)*hp(1).field(EZ,z))
		    /invomep0*(-1.0)*epseff.eps(z);
		v = vp.field(HY, x);
		c = CCI*invomep0/str.eps(x, z);
		return c*h*v;
		break;
	case EY: return CC0;
		break;
	case EZ: 
		h = amp(0)*hp(0).field(HY,z)+amp(1)*hp(1).field(HY,z);
		v = vp.field(EZ, x)/invomep0*(-1.0)*vp.wg.eps(x);
		c = CCmI*invomep0/str.eps(x, z);
		return c*h*v;
		break;
	case HX: return CC0;
		break;
	case HY: 
		h = amp(0)*hp(0).field(HY,z)+amp(1)*hp(1).field(HY,z);
		return h*vp.field(HY, x);
		break;
	case HZ: return CC0;
		break;
	default: eimerror("field: cp not implemented"); 
		break;
	}
	return CC0;
}

/* values on a rectangular mesh */
Cmatrix EimField::field(Fcomp cp, Interval wx, int npx, 
                                   Interval wz, int npz) const
{
	Cmatrix f(npx, npz);
	double dx = (wx.x1-wx.x0)/((double)(npx-1));
	double dz = (wz.x1-wz.x0)/((double)(npz-1));
	double x, z;
	z = wz.x0;
	for(int j=0; j<=npz-1; ++j)
	{
		x = wx.x0;
		for(int i=0; i<=npx-1; ++i) 
		{ 
			f(i, j) = field(cp, x, z); 
			x += dx;
		}
		z += dz;
	}
	return f;
}

/* multiply the entire field by a phase factor exp(i phi),
   the internal representaion is modified accordingly */
void EimField::adjustphase(double phi)
{
	Complex a;
	a = expi(phi);
	amp.multeq(a);
	return;
}

/* adjust the global phase, such that a plot(cp, ORG) of the physical
   basic field component (cp = EY for TE, cp = HY for TM) exhibits the
   maximum field at position (x, z) */
void EimField::adjustphase(double x, double z)
{
	Fcomp cp = principalcomp(pol);
	adjustphase(-field(cp, x, z).arg());
	return;
}

/* adjust the global phase, such that a plot(cp, ORG) of the physical
   basic field component (cp = EY for TE, cp = HY for TM) exhibits an 
   overall field maximum (-> time snapshot of an extremal state  
   in case of a configuration with partially standing waves)   
   ... determine the maximum on a regular mesh of numx * mumz points
       on the window given by xbeg, xend, zbeg, zend */ 
void EimField::adjustphase(double xbeg, double xend, double zbeg, double zend,
                            int numx, int numz)
{
	Fcomp cp = principalcomp(pol);
	double x, z;
	double xs, zs;
	xs = (xend-xbeg)/((double) numx);
	zs = (zend-zbeg)/((double) numz);
	x = xbeg;
	double xm = 0.0, zm = 0.0, maxf = 0.0;
	double f;
	for(int xi=0; xi <=numx; ++xi)
	{
		z = zbeg;
		for(int zi=0; zi <=numz; ++zi)
		{
			f = field(cp, x, z).sqabs();
			if(f > maxf) { xm=x; zm=z; maxf=f; }
			z += zs;
		}
		x += xs;
	}
	adjustphase(-field(cp, xm, zm).arg());
	return;
}

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

/* mode interference plots, image of component cp
   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 EimField::plot(Fcomp cp, Afo foa,
                     double xbeg, double xend, double zbeg, double zend, 
                     int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "pl_____I.m";
	name[2] = afochr(foa);
        name[3] = fldchr(cp);
        name[4] = cpchr(cp);      
	name[5] = ext0;
	name[6] = ext1;
	fprintf(stderr, "plot([%g, %g] x [%g, %g]) >> %s\n", 
	                xbeg, xend, zbeg, zend, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix cfld = field(cp, Interval(xbeg, xend), npx, 
	                         Interval(zbeg, zend), npz);
	Dmatrix fld = realize(cfld, foa); 
	mlout_fld(dat, npx, npz, cp, fld);
	name[8] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage(cp, foa, name, desc, dat); 
	if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
	else                         mlout_magcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* export full field data into a viewer m-file
   xbeg, xend, zbeg, zend: x resp. z extensions of the plot window
   npx, npz: number of plot points in the x and z directions
   ext0, ext1: filename id characters */
void EimField::viewer(double xbeg, double xend, double zbeg, double zend,
                       int npx, int npz, char ext0, char ext1) const
{
	FILE *dat;
	char name[13] = "fld__A.m";
	name[3] = ext0;
	name[4] = ext1;
	fprintf(stderr, 
	   "viewer([%g (%d) %g] x [%g (%d) %g]) >> %s\n", 
	   xbeg, npx, xend, zbeg, npz, zend, name);
	dat = fopen(name, "w+");
	name[6] = 0; 
	mlout_viewertop(dat, name, pol, lambda);
	mlout_gengeoxz(dat, str, xbeg, xend, zbeg, zend);
	mlout_meshxz(dat, xbeg, xend, npx, zbeg, zend, npz);
	Cmatrix f;
	Fcomp cp;
	if(pol == TE)
	{
		cp = EX;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EY;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EZ;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HX;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HY;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HZ;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
	}
	else 	
	{
		cp = EX;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = EY;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = EZ;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HX;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
		cp = HY;
		f = field(cp, Interval(xbeg, xend), npx, 
			      Interval(zbeg, zend), npz); 
		mlout_fld(dat, npx, npz, cp, f.re());
		mlout_fldtore(dat, cp);
		mlout_fld(dat, npx, npz, cp, f.im());
		mlout_fldtoim(dat, cp);
		cp = HZ;
		mlout_0fld(dat, cp);
		mlout_fldtore(dat, cp);
		mlout_0fld(dat, cp);
		mlout_fldtoim(dat, cp);
	}

	Dmatrix n(npx, npz);
	double dx = (xend-xbeg)/((double)(npx-1));
	double dz = (zend-zbeg)/((double)(npz-1));
	double x, z;
	z = zbeg;
	for(int j=0; j<=npz-1; ++j)
	{
		x = xbeg;
		for(int i=0; i<=npx-1; ++i) { n(i, j) = str.n(x, z); x += dx; }
		z += dz;
	}
	mlout_n(dat, npx, npz, n);
	mlout_fldviewer(dat, name);
	fclose(dat);
	return;
}