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

/*
 * matlvis.cpp
 * Visualization, helper routines, output to MATLAB m-files
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"gengwed.h"
#include"structure.h"
#include"matlvis.h"

/* some global settings */

/* output uses colour */
int Mlop_Colour=YES;
/* output files include print commands */
int Mlop_Print=NO;
/* add contour lines to image plots */
int Mlop_Contour=NO;

/* linear black & blue colormap */
void mlout_cmlinblackblue(FILE *dat)
{
	fprintf(dat, "%% Choose colormap\n");
	fprintf(dat, "mi = 32; mlv = 0.2;\n");
	fprintf(dat, "blue = zeros(64, 3);\n");
	fprintf(dat, "blue(1:mi, 3) = linspace(0, 1,   mi)';\n");
	fprintf(dat, "blue(1:mi, 2) = linspace(0, mlv, mi)';\n");
	fprintf(dat, "blue(1:mi, 1) = linspace(0, mlv, mi)';\n");
	fprintf(dat, "blue(mi+1:64, 3) = 1;\n");
	fprintf(dat, "blue(mi+1:64, 2) = linspace(mlv, 1, 64-mi)';\n");
	fprintf(dat, "blue(mi+1:64, 1) = linspace(mlv, 1, 64-mi)'; \n");
	fprintf(dat, "colormap(blue);\n");
	fprintf(dat, "\n");
	return;
}

/* linear blue colormap */
void mlout_cmlinblue(FILE *dat)
{
	fprintf(dat, "%% Choose colormap\n");
	fprintf(dat, "blue = zeros(64, 3);\n");
	fprintf(dat, "blue(:,1) = linspace(0,1,64)';\n");
	fprintf(dat, "blue(:,2) = blue(:,1);\n");
	fprintf(dat, "blue(:,3) = 1;\n");
	fprintf(dat, "colormap(blue);\n");
	fprintf(dat, "\n");
	return;
}

/* linear black colormap */
void mlout_cmlinblack(FILE *dat)
{
	fprintf(dat, "colormap(gray);\n");
	fprintf(dat, "\n");
	return;
}

/* black & blue colormap, shows low levels in more detail */
void mlout_cmlowlevblackblue(FILE *dat)
{
	fprintf(dat, "%% Choose colormap\n");
	fprintf(dat, "blue = zeros(256, 3);\n");
	fprintf(dat, 
	    "cpd = 0.01; cphb= 0.75; cpho= 0.25; cpa = atan(cpd*256);\n");
	fprintf(dat, "for i = 0:255\n");
	fprintf(dat, "   t = atan(cpd*i)/cpa;\n");
	fprintf(dat, "   blue(i+1, 3) = (2-4*cphb)*t*t+(4*cphb-1)*t;\n");
	fprintf(dat, "   blue(i+1, 1) = (2-4*cpho)*t*t+(4*cpho-1)*t;\n");
	fprintf(dat, "   blue(i+1, 2) = blue(i+1, 1);\n");
	fprintf(dat, "end\n");
	fprintf(dat, "colormap(blue);\n");
	fprintf(dat, "\n");
	return;
}

/* black & blue colormap, magnified center */
void mlout_cmmagblackblue(FILE *dat)
{
	fprintf(dat, "%% Choose colormap\n");
	fprintf(dat, "blue = zeros(256, 3);\n");
	fprintf(dat, 
	    "cpd = 0.02; cphb= 0.75; cpho= 0.25; cpa = atan(256*cpd/2);\n");
	fprintf(dat, "for i = 0:255\n");
	fprintf(dat, "   t = 0.5*(atan(cpd*(i-256/2))/cpa+1);\n");
	fprintf(dat, "   blue(i+1, 3) = (2-4*cphb)*t*t+(4*cphb-1)*t;\n");
	fprintf(dat, "   blue(i+1, 1) = (2-4*cpho)*t*t+(4*cpho-1)*t;\n");
	fprintf(dat, "   blue(i+1, 2) = blue(i+1, 1);\n");
	fprintf(dat, "end\n");
	fprintf(dat, "colormap(blue);\n");
	fprintf(dat, "\n");
	return;
}

/* linear colormap */
void mlout_lincolormap(FILE *dat)
{
	if(Mlop_Colour==YES) mlout_cmlinblackblue(dat);
	else                 mlout_cmlinblack(dat);
	return;
}

/* colormap for positive values, shows low levels in more detail */
void mlout_lowlevcolormap(FILE *dat)
{
	if(Mlop_Colour==YES) mlout_cmlowlevblackblue(dat);
	else                 mlout_cmlinblack(dat);
	return;
}

/* colormap, magnified range around the center (only for Colour=YES) */
void mlout_magcolormap(FILE *dat)
{
	if(Mlop_Colour==YES) mlout_cmmagblackblue(dat);
	else                 mlout_cmlinblack(dat);
	return;
}

/* sketch of permittivity boundaries: colours, linewidth */
void mlout_sketchinit(FILE *dat)
{
	fprintf(dat, "%% Add sketch of permittivity boundaries\n");
	fprintf(dat, "hold on;\n");
	fprintf(dat, "lw = 1.0;\n");
	if(Mlop_Colour==YES) fprintf(dat, "ls = 'w-';\n");
	else                 fprintf(dat, "ls = 'w-';\n");
	return;
}

/* sketch of permittivity boundaries: outline the geometry */
void mlout_geooutline(FILE *dat)
{
	fprintf(dat, "for i=1:nhd\n");
	fprintf(dat, "   plot([hdbeg(i), hdend(i)], [hdpos(i), hdpos(i)], ls, 'LineWidth', lw);\n");
	fprintf(dat, "end\n");
	fprintf(dat, "for i=1:nvd\n");
	fprintf(dat, "   plot([vdpos(i), vdpos(i)], [vdbeg(i), vdend(i)], ls, 'LineWidth', lw);\n");
	fprintf(dat, "end\n");
	fprintf(dat, "ra = [0:2*pi/500:2*pi];\n");
	fprintf(dat, "for i=1:ncd\n");
	fprintf(dat, "   plot(cos(ra)*cdrad(i)+cdhp(i), sin(ra)*cdrad(i)+cdvp(i), ...\n");
	fprintf(dat, "        ls, 'LineWidth', lw);\n");
	fprintf(dat, "end\n");
	fprintf(dat, "\n");
	return;
}

/* output window, x-z-plane */
void mlout_windowxz(FILE *dat, double xbeg, double xend,
                               double zbeg, double zend)
{
	mlout_window(dat, 'x', xbeg, xend, 'z', zbeg, zend);
	return;
}

/* output window, generic */
void mlout_window(FILE *dat, 
                  char xchr, double xbeg, double xend,
                  char zchr, double zbeg, double zend)
{
	fprintf(dat, "%% Output window\n");
	fprintf(dat, "%cbeg = %g;\n", xchr, xbeg);
	fprintf(dat, "%cend = %g;\n", xchr, xend);
	fprintf(dat, "%cbeg = %g;\n", zchr, zbeg);
	fprintf(dat, "%cend = %g;\n", zchr, zend);
	fprintf(dat, "\n");
	return;
}

/* output mesh, x-z-plane */
void mlout_meshxz(FILE *dat, double xbeg, double xend, int npx,
                             double zbeg, double zend, int npz)
{
	mlout_mesh(dat, 'x', xbeg, xend, npx, 'z', zbeg, zend, npz);
	return;
}

/* output mesh, generic */
void mlout_mesh(FILE *dat, 
                char xchr, double xbeg, double xend, int npx,
                char zchr, double zbeg, double zend, int npz)
{
	mlout_window(dat, xchr, xbeg, xend, zchr, zbeg, zend);
	fprintf(dat, "%% Mesh information\n");
	fprintf(dat, "num%c = %d;\n", xchr, npx);
	if(xbeg==xend)
		fprintf(dat, "%c = %g*ones(num%c, 1);\n", xchr, xbeg, xchr);
	else
		fprintf(dat, "%c = linspace(%cbeg, %cend, num%c);\n", xchr, xchr, xchr, xchr);
	fprintf(dat, "num%c = %d;\n", zchr, npz);
	if(zbeg==zend)
		fprintf(dat, "%c = %g*ones(num%c, 1);\n", zchr, zbeg, zchr);
	else
		fprintf(dat, "%c = linspace(%cbeg, %cend, num%c);\n", zchr, zchr, zchr, zchr);
	fprintf(dat, "\n");
	return;
}

/* previously stored field is meant as real part of a field component */
void mlout_fldtore(FILE *dat, Fcomp cp)
{
fprintf(dat, "%% Meant as the real part\n");
fprintf(dat, "%c%c_re = %c%c;\n", fldchr(cp),cpchr(cp),fldchr(cp),cpchr(cp)); 
fprintf(dat, "\n");
}

/* previously stored field is meant as imaginary part of a field component */
void mlout_fldtoim(FILE *dat, Fcomp cp)
{
fprintf(dat, "%% Meant as the imaginary part\n");
fprintf(dat, "%c%c_im = %c%c;\n", fldchr(cp),cpchr(cp),fldchr(cp),cpchr(cp)); 
fprintf(dat, "\n");
}

/* general structure: geometry for interference visualization */
void mlout_gengeoxz(FILE *dat, SegWgStruct s, 
                    double xbeg, double xend, double zbeg, double zend) 
{
	fprintf(dat, "%% Geometry\n");
	fprintf(dat, "\n");

	Dvector xpos;
	Dvector beg;
	Dvector end;
	Dvector pos;
	int x1, x2, p;
	double x1p, x2p, z0p, z1p, r; 
	int nd;

	fprintf(dat, "%% Vertical discontinuities\n");
	for(int j=0; j<=s.nz; ++j)
	{
		x1 = 0; x2 = 0; xpos.strip();
		while(x1 <= s(j).nx || x2 <= s(j+1).nx)
		{
			if(x1<=s(j).nx) x1p = s(j).hx(x1); 
			else x1p=1.0e+300;
			if(x2<=s(j+1).nx) x2p = s(j+1).hx(x2); 
			else x2p=1.0e+300;
			if(x1p < x2p)
			{
				xpos.append(x1p);
				++x1;
			}
			else
			{
				xpos.append(x2p);
				++x2;
			}
		}
		p = xpos.nel;
		for(int i=-1; i<=p-1; ++i)
		{
			if(i==-1) x1p = xbeg; 
			else      x1p = xpos(i);
			if(i==p-1) x2p = xend;
			else       x2p = xpos(i+1);
			if(fabs(x1p-x2p) >= 1.0e-10)
			{
				r = 0.5*(x1p+x2p);
				if(fabs(s(j).eps(r)-s(j+1).eps(r))
				   >1.0e-10)
				{
					beg.append(x1p);
					end.append(x2p);
					pos.append(s.hz(j));
				} 
			}
		}
	}
	nd = pos.nel;
	fprintf(dat, "nvd = %d;\n", nd);
	fprintf(dat, "vdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	pos.strip();
	beg.strip();
	end.strip();

	fprintf(dat, "%% Horizontal discontinuities\n");
	for(int j=0; j<=s.nz+1; ++j)
	{
		if(j == 0) z0p = zbeg; 
		else       z0p = s.hz(j-1);
		if(j == s.nz+1) z1p = zend; 
		else       	z1p = s.hz(j);
		for(int i=0; i<=s(j).nx; ++i)
		{
			if(fabs(s(j).eps(i)-s(j).eps(i+1))>1.0e-10)
			{
				pos.append(s(j).hx(i));
				beg.append(z0p);
				end.append(z1p);
			}
		}
	}
	nd = pos.nel;
	fprintf(dat, "nhd = %d;\n", nd);
	fprintf(dat, "hdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	fprintf(dat, "%% Circular discontinuities\n");
	fprintf(dat, "ncd = %d;\n", 0);
	fprintf(dat, "cdrad = [ ");
	fprintf(dat, "];\n");
	fprintf(dat, "cdvp  = [ ");
	fprintf(dat, "];\n");
	fprintf(dat, "cdhp  = [ ");
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	return;
}


/* overlay structure: geometry for interference visualization */
void mlout_gengeoxz(FILE *dat, OvlStruct o, 
                    double xbeg, double xend, double zbeg, double zend) 
{
	fprintf(dat, "%% Geometry\n");
	fprintf(dat, "\n");

	Dvector beg;
	Dvector end;
	Dvector pos;
	Dvector rad;
	Dvector x0;
	Dvector z0;
	int nd;

	fprintf(dat, "%% Vertical discontinuities\n");
	for(int j=0; j<=o.no-1; ++j)
	{
		switch(o(j).type)
		{
		case VLAY:
			pos.append(o(j).ubd);
			beg.append(xbeg);
			end.append(xend);
			pos.append(o(j).lbd);
			beg.append(xbeg);
			end.append(xend);
			break;	
		case HLAY: case RING: case DISK: case RECT:
			break;	
		}
	}
	nd = pos.nel;
	fprintf(dat, "nvd = %d;\n", nd);
	fprintf(dat, "vdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	pos.strip();
	beg.strip();
	end.strip();

	fprintf(dat, "%% Horizontal discontinuities\n");
	for(int j=0; j<=o.no-1; ++j)
	{
		switch(o(j).type)
		{
		case HLAY:
			pos.append(o(j).ubd);
			beg.append(zbeg);
			end.append(zend);
			pos.append(o(j).lbd);
			beg.append(zbeg);
			end.append(zend);
			break;	
		case VLAY: case RING: case DISK: case RECT:
			break;	
		}
	}
	nd = pos.nel;
	fprintf(dat, "nhd = %d;\n", nd);
	fprintf(dat, "hdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	fprintf(dat, "%% Circular discontinuities\n");
	for(int j=0; j<=o.no-1; ++j)
	{
		switch(o(j).type)
		{
		case RING:
			rad.append(o(j).ubd);
			x0.append(o(j).xo);
			z0.append(o(j).zo);
			rad.append(o(j).lbd);
			x0.append(o(j).xo);
			z0.append(o(j).zo);
			break;	
		case DISK:
			rad.append(o(j).ubd);
			x0.append(o(j).xo);
			z0.append(o(j).zo);
			break;	
		case HLAY: case VLAY: case RECT:
			break;	
		}
	}
	nd = rad.nel;
	fprintf(dat, "ncd = %d;\n", nd);
	fprintf(dat, "cdrad = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", rad(i));
	fprintf(dat, "];\n");
	fprintf(dat, "cdvp  = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", x0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "cdhp  = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", z0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	return;
}
void mlout_gengeoxy(FILE *dat, OvlStruct o, Rect disp) 
{
	fprintf(dat, "%% Geometry\n");
	fprintf(dat, "\n");

	Dvector beg;
	Dvector end;
	Dvector pos;
	Dvector rad;
	Dvector x0;
	Dvector y0;
	int nd;

	fprintf(dat, "%% Vertical discontinuities\n");
	for(int j=0; j<=o.no-1; ++j)
	{
		switch(o(j).type)
		{
		case VLAY:
			pos.append(o(j).ubd);
			beg.append(disp.x0);
			end.append(disp.x1);
			pos.append(o(j).lbd);
			beg.append(disp.x0);
			end.append(disp.x1);
			break;	
		case HLAY: case RING: case DISK: case RECT:
			break;	
		}
	}
	nd = pos.nel;
	fprintf(dat, "nvd = %d;\n", nd);
	fprintf(dat, "vdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	pos.strip();
	beg.strip();
	end.strip();

	fprintf(dat, "%% Horizontal discontinuities\n");
	for(int j=0; j<=o.no-1; ++j)
	{
		switch(o(j).type)
		{
		case HLAY:
			pos.append(o(j).ubd);
			beg.append(disp.y0);
			end.append(disp.y1);
			pos.append(o(j).lbd);
			beg.append(disp.y0);
			end.append(disp.y1);
			break;	
		case VLAY: case RING: case DISK: case RECT:
			break;	
		}
	}
	nd = pos.nel;
	fprintf(dat, "nhd = %d;\n", nd);
	fprintf(dat, "hdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	fprintf(dat, "%% Circular discontinuities\n");
	for(int j=0; j<=o.no-1; ++j)
	{
		switch(o(j).type)
		{
		case RING:
			rad.append(o(j).ubd);
			x0.append(o(j).xo);
			y0.append(o(j).zo);
			rad.append(o(j).lbd);
			x0.append(o(j).xo);
			y0.append(o(j).zo);
			break;	
		case DISK:
			rad.append(o(j).ubd);
			x0.append(o(j).xo);
			y0.append(o(j).zo);
			break;	
		case HLAY: case VLAY: case RECT:
			break;	
		}
	}
	nd = rad.nel;
	fprintf(dat, "ncd = %d;\n", nd);
	fprintf(dat, "cdrad = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", rad(i));
	fprintf(dat, "];\n");
	fprintf(dat, "cdvp  = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", x0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "cdhp  = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", y0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	return;
}


/* general structure: refractive index profile patches, data */
void mlout_gengeorefindpatches(FILE *dat, SegWgStruct s) 
{
        fprintf(dat, "%% Refractive index patches\n");
	int sp = s.special();
	Circuit c = s.circuit();
	double faraway = (c.hx(c.nx)-c.hx(0)+c.hz(c.nz)-c.hz(0))*10.0;
	int nump = (c.nx+2)*(c.nz+2);
	Dvector vx0(nump), vx1(nump), vz0(nump), vz1(nump), vn(nump);
	nump = 0;
	double x0, x1, z0, z1;
	for(int xi=0; xi<=c.nx+1; ++xi)
	{
		if(xi == 0) x0 = -faraway;
		else        x0 = c.hx(xi-1);
		if(xi == c.nx+1) x1 = faraway;
		else             x1 = c.hx(xi);
		for(int zi=0; zi<=c.nz+1; ++zi)
		{
			if(zi == 0) z0 = -faraway;
			else        z0 = c.hz(zi-1);
			if(zi == c.nz+1) z1 = faraway;
			else             z1 = c.hz(zi);
			vx0(nump) = x0;
			vx1(nump) = x1;
			vz0(nump) = z0;
			vz1(nump) = z1;
			if(sp) vn(nump) = sgn_sqrt(c.n(xi, zi)); 
			else   vn(nump) = c.n(xi, zi); 
			++nump;
		}
	}
	fprintf(dat, "numnp = %d;\n", nump);
	fprintf(dat, "npx0 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vx0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npx1 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vx1(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npz0 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vz0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npz1 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vz1(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npn = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vn(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
}

/* waveguide geometry and field extrema */
void mlout_geo(FILE *dat, Waveguide wg, double minf, double maxf)
{
        int i;
        fprintf(dat, "%% Waveguide geometry\n");
        fprintf(dat, "nlX = %d;\n", wg.nx);
        fprintf(dat, "bdX = [");
        for(i=0; i<=wg.nx; ++i) fprintf(dat, "%g ", wg.hx(i));
        fprintf(dat, "];\n");
        fprintf(dat, "\n");
        fprintf(dat, "%% Refractive indices\n");
        fprintf(dat, "n = [");
	if(wg.special)
	{ 
        	for(i=0; i<=wg.nx+1; ++i) 
			fprintf(dat, "%g ", sgn_sqrt(wg.n(i)));
	}
	else
	{ 
        	for(i=0; i<=wg.nx+1; ++i) 
			fprintf(dat, "%g ", wg.n(i));
	}
        fprintf(dat, "];\n");
        fprintf(dat, "\n");
        fprintf(dat, "%% Bounds on field values\n");
        fprintf(dat, "minf = %g;\n", minf);
        fprintf(dat, "maxf = %g;\n", maxf);
        fprintf(dat, "\n");
        return;
}

/* add contour lines to an image plot */
void mlout_contours(Fcomp cp, FILE *dat)
{
fprintf(dat, "%% Add contour lines\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "numct = 10;\n");
fprintf(dat, "if minf > -(maxf-minf)/numct\n");
fprintf(dat, "   cl0 = minf;\n");
fprintf(dat, "   cl1 = maxf;\n");
fprintf(dat, "else\n");
fprintf(dat, "   cl1 = max([abs(maxf), abs(minf)]);\n");
fprintf(dat, "   cl0 = -cl1;\n");
fprintf(dat, "end\n");
fprintf(dat, "lw = 0.5;\n");
fprintf(dat, "ls = '-';\n");
fprintf(dat, "lc = 'b';\n");
fprintf(dat, "for ci=1:numct;\n");
fprintf(dat, "   lv = cl0+ci*(cl1-cl0)/(numct+1);\n");
fprintf(dat, "   lim = 0;\n");
	if(Mlop_Contour == YES)
	{
fprintf(dat, "   cm = contourc(z, x, %c%c, [lv, lv]); lim = size(cm, 2);\n",
		fldchr(cp), cpchr(cp));
	}
	else
	{
fprintf(dat, "%%   cm = contourc(z, x, %c%c, [lv, lv]); lim = size(cm, 2);\n",
		fldchr(cp), cpchr(cp));
	}
fprintf(dat, "   i = 1;\n");
fprintf(dat, "   while(i < lim)\n");
fprintf(dat, "      np = cm(2, i);\n");
fprintf(dat, "      pz = cm(1, i+1:i+np);\n");
fprintf(dat, "      px = cm(2, i+1:i+np);\n");
fprintf(dat, "      line(pz, px, 'Color', lc, 'LineStyle', ls, 'LineWidth', lw);\n");
fprintf(dat, "      i = i+np+1;\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* add contour lines to an image movie */
void mlout_contours(FILE *dat)
{
fprintf(dat, "%% Add contour lines\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "numct = 10;\n");
fprintf(dat, "lw = 0.5;\n");
fprintf(dat, "ls = '-';\n");
fprintf(dat, "lc = 'b';\n");
fprintf(dat, "for ci=1:numct;\n");
fprintf(dat, "   lv = -maxf+ci*2*maxf/(numct+1);\n");
fprintf(dat, "   lim = 0;\n");
	if(Mlop_Contour == YES)
	{
fprintf(dat, "   cm = contourc(z, x, fld, [lv, lv]); lim = size(cm, 2);\n");
	}
	else
	{
fprintf(dat, "%%   cm = contourc(z, x, fld, [lv, lv]); lim = size(cm, 2);\n");
	}
fprintf(dat, "   i = 1;\n");
fprintf(dat, "   while(i < lim)\n");
fprintf(dat, "      np = cm(2, i);\n");
fprintf(dat, "      pz = cm(1, i+1:i+np);\n");
fprintf(dat, "      px = cm(2, i+1:i+np);\n");
fprintf(dat, "      line(pz, px, 'Color', lc, 'LineStyle', ls, 'LineWidth', lw);\n");
fprintf(dat, "      i = i+np+1;\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* add contour lines to an image plot, generic */
void mlout_contours(char xchr, char zchr, Fcomp cp, FILE *dat)
{
fprintf(dat, "%% Add contour lines\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "numct = 10;\n");
fprintf(dat, "if minf > -(maxf-minf)/numct\n");
fprintf(dat, "   cl0 = minf;\n");
fprintf(dat, "   cl1 = maxf;\n");
fprintf(dat, "else\n");
fprintf(dat, "   cl1 = max([abs(maxf), abs(minf)]);\n");
fprintf(dat, "   cl0 = -cl1;\n");
fprintf(dat, "end\n");
fprintf(dat, "lw = 0.5;\n");
fprintf(dat, "ls = '-';\n");
fprintf(dat, "lc = 'b';\n");
fprintf(dat, "for ci=1:numct;\n");
fprintf(dat, "   lv = cl0+ci*(cl1-cl0)/(numct+1);\n");
fprintf(dat, "   lim = 0;\n");
	if(Mlop_Contour == YES)
	{
fprintf(dat, "   cm = contourc(%c, %c, %c%c, [lv, lv]); lim = size(cm, 2);\n",
		zchr, xchr, fldchr(cp), cpchr(cp));
	}
	else
	{
fprintf(dat, "%%   cm = contourc(%c, %c, %c%c, [lv, lv]); lim = size(cm, 2);\n",
		zchr, xchr, fldchr(cp), cpchr(cp));
	}
fprintf(dat, "   i = 1;\n");
fprintf(dat, "   while(i < lim)\n");
fprintf(dat, "      np = cm(2, i);\n");
fprintf(dat, "      pz = cm(1, i+1:i+np);\n");
fprintf(dat, "      px = cm(2, i+1:i+np);\n");
fprintf(dat, "      line(pz, px, 'Color', lc, 'LineStyle', ls, 'LineWidth', lw);\n");
fprintf(dat, "      i = i+np+1;\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* add contour lines to an image movie, generic */
void mlout_contours(char xchr, char zchr, FILE *dat)
{
fprintf(dat, "%% Add contour lines\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "numct = 10;\n");
fprintf(dat, "lw = 0.5;\n");
fprintf(dat, "ls = '-';\n");
fprintf(dat, "lc = 'b';\n");
fprintf(dat, "for ci=1:numct;\n");
fprintf(dat, "   lv = -maxf+ci*2*maxf/(numct+1);\n");
fprintf(dat, "   lim = 0;\n");
	if(Mlop_Contour == YES)
	{
fprintf(dat, "   cm = contourc(%c, %c, fld, [lv, lv]); lim = size(cm, 2);\n", zchr, xchr);
	}
	else
	{
fprintf(dat, "%%   cm = contourc(%c, %c, fld, [lv, lv]); lim = size(cm, 2);\n", zchr, xchr);
	}
fprintf(dat, "   i = 1;\n");
fprintf(dat, "   while(i < lim)\n");
fprintf(dat, "      np = cm(2, i);\n");
fprintf(dat, "      pz = cm(1, i+1:i+np);\n");
fprintf(dat, "      px = cm(2, i+1:i+np);\n");
fprintf(dat, "      line(pz, px, 'Color', lc, 'LineStyle', ls, 'LineWidth', lw);\n");
fprintf(dat, "      i = i+np+1;\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* mode section plot geometry  */
void mlout_Lsecgeo(FILE *dat, double x0, double x1,
                   int nbd, Dvector bd, Dvector ri)
{
        int i;
fprintf(dat, "%% Section location\n");
fprintf(dat, "x0 = %g;\n", x0);
fprintf(dat, "x1 = %g;\n", x1);
fprintf(dat, "\n");
fprintf(dat, "%% Crossed boundaries \n");
fprintf(dat, "nbd = %d;\n", nbd);
fprintf(dat, "bd = [");
        for(i=0; i<=nbd-1; ++i)
        {
fprintf(dat, "%g ", bd(i));
        }
fprintf(dat, "];\n");
fprintf(dat, "ri = [");
        for(i=0; i<=nbd; ++i)
        {
fprintf(dat, "%g ", ri(i));
        }
fprintf(dat, "];\n");
fprintf(dat, "\n");
        return;
}

/* section of a 2D field, begin + end points */
void mlout_sec(FILE *dat, double x0, double z0, double x1, double z1, 
               int np, Fcomp cp, int nc, Dvector fld)
{
	fprintf(dat, "%c%c%d = [", fldchr(cp), cpchr(cp), nc);
	for(int j=0; j<=np-1; ++j) fprintf(dat, "%g ", fld(j));
	fprintf(dat, "]; ");
	fprintf(dat, "xp%d = linspace(%g, %g, %d); ", nc, x0, x1, np);
	fprintf(dat, "zp%d = linspace(%g, %g, %d);\n", nc, z0, z1, np);
	return;
}

/* section of a complex 2D field, begin + end points */
void mlout_sec(FILE *dat, double x0, double z0, double x1, double z1, 
               int np, Fcomp cp, int nc, Cvector fld)
{
	Dvector v;
	v = fld.re();
	fprintf(dat, "%c%c%d_re = [", fldchr(cp), cpchr(cp), nc);
	for(int j=0; j<=np-1; ++j) fprintf(dat, "%g ", v(j));
	fprintf(dat, "]; ");
	v = fld.im();
	fprintf(dat, "%c%c%d_im = [", fldchr(cp), cpchr(cp), nc);
	for(int j=0; j<=np-1; ++j) fprintf(dat, "%g ", v(j));
	fprintf(dat, "]; ");
	fprintf(dat, "xp%d = linspace(%g, %g, %d); ", nc, x0, x1, np);
	fprintf(dat, "zp%d = linspace(%g, %g, %d);\n", nc, z0, z1, np);
	return;
}

/* for the fancy plots: transfer data related to the geometry annotations */
void mlout_bddata(FILE *dat, Fcomp cp, int nc, Ivector si, Ivector np, 
                  Dvector f, Dvector x0, Dvector z0, Dvector x1, Dvector z1)
{
	fprintf(dat, "%% Interface section data\n");
	fprintf(dat, "nc = %d;\n", nc);
	fprintf(dat, "bdsi = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%d ", si(j)+1);
	fprintf(dat, "];\n");
	fprintf(dat, "bdnp = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%d ", np(j));
	fprintf(dat, "];\n");

	fprintf(dat, "bdf%c%c = [", fldchr(cp), cpchr(cp));
	for(int s=0; s<=nc-1; ++s)
	{
		int p0 = si(s);
		for(int p=0; p<=np(s)-1; ++p) fprintf(dat, "%g ", f(p0+p));
	}
	fprintf(dat, "];\n");

	fprintf(dat, "bdx0 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", x0(j));
	fprintf(dat, "];\n");
	fprintf(dat, "bdz0 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", z0(j));
	fprintf(dat, "];\n");
	fprintf(dat, "bdx1 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", x1(j));
	fprintf(dat, "];\n");
	fprintf(dat, "bdz1 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", z1(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	return;
	
}
void mlout_bddata(FILE *dat, Fcomp cp, int nc, Ivector si, Ivector np, 
                  Cvector f, Dvector x0, Dvector z0, Dvector x1, Dvector z1)
{
	fprintf(dat, "%% Interface section data\n");
	fprintf(dat, "nc = %d;\n", nc);
	fprintf(dat, "bdsi = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%d ", si(j)+1);
	fprintf(dat, "];\n");
	fprintf(dat, "bdnp = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%d ", np(j));
	fprintf(dat, "];\n");

	fprintf(dat, "bdf%c%c_re = [", fldchr(cp), cpchr(cp));
	for(int s=0; s<=nc-1; ++s)
	{
		int p0 = si(s);
		for(int p=0; p<=np(s)-1; ++p) fprintf(dat, "%g ", f(p0+p).re);
	}
	fprintf(dat, "];\n");

	fprintf(dat, "bdf%c%c_im = [", fldchr(cp), cpchr(cp));
	for(int s=0; s<=nc-1; ++s)
	{
		int p0 = si(s);
		for(int p=0; p<=np(s)-1; ++p) fprintf(dat, "%g ", f(p0+p).im);
	}
	fprintf(dat, "];\n");

	fprintf(dat, "bdx0 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", x0(j));
	fprintf(dat, "];\n");
	fprintf(dat, "bdz0 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", z0(j));
	fprintf(dat, "];\n");
	fprintf(dat, "bdx1 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", x1(j));
	fprintf(dat, "];\n");
	fprintf(dat, "bdz1 = [");
	for(int j=0; j<=nc-1; ++j) fprintf(dat, "%g ", z1(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	return;
}

/* 1D field section, including position vector */
void mlout_sec(FILE *dat,
               Fcomp cp, char c2, char c3, char c4, char c5,
               int np, Dvector fld, Dvector pos)
{
        int j;

        fprintf(dat, "%c%c%c%c%c%c = [",
                     fldchr(cp), cpchr(cp), c2, c3, c4, c5);
        for(j=0; j<=np-1; ++j) fprintf(dat, "%g ", fld(j));
        fprintf(dat, "];\n");
        fprintf(dat, "pos%c%c%c%c = [", c2, c3, c4, c5);
        for(j=0; j<=np-1; ++j) fprintf(dat, "%g ", pos(j));
        fprintf(dat, "];\n");
        return;
}

/* mode section plot, plot commands */
void mlout_Lsecplot(const char *name, FILE *dat, 
                    Fcomp cp, Afo foa, int nbd, int nsec)
{
	int i;
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "axis off;\n");
fprintf(dat, "alw = 1.2;\n");
fprintf(dat, "lw = 1.6;\n");
fprintf(dat, "fs = 16;\n");
fprintf(dat, "lfs = 16;\n");
fprintf(dat, "fnm = 'Helvetica';\n");
fprintf(dat, "\n");

fprintf(dat, "y0 = minf-(maxf-minf)/20.0;\n");
fprintf(dat, "y1 = maxf+(maxf-minf)/20.0;\n");
fprintf(dat, "%% y0 = minf;\n");
fprintf(dat, "%% y1 = maxf;\n");
fprintf(dat, "\n");

fprintf(dat, "h = axes('Position', [0.15 0.15 0.80 0.65], ...\n");
fprintf(dat, "       'XLim', [x0 x1], ...\n");
fprintf(dat, "       'YLim', [y0 y1], ...\n");
fprintf(dat, "       'FontSize', fs, ...\n");
fprintf(dat, "       'FontName', fnm, ...\n");
fprintf(dat, "       'LineWidth', alw, ...\n");
fprintf(dat, "       'TickDir', 'out');\n");

fprintf(dat, "axes(h);\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");

fprintf(dat, "xlabel(['x [' char(956) 'm]'], 'FontSize', lfs, 'FontName', fnm);\n");
switch(foa)
{
case MOD:
	fprintf(dat, "ylabel('|%c_%c|', 'FontSize', lfs, 'Fontname', fnm);\n",
                fldchr(cp), cpchr(cp));
	break;
case SQR: 
	fprintf(dat, "ylabel('|%c_%c|^2', 'FontSize', lfs, 'Fontname', fnm);\n",
                fldchr(cp), cpchr(cp));
	break;
case ORG:
	fprintf(dat, "ylabel('%c_%c', 'FontSize', lfs, 'Fontname', fnm);\n",
                fldchr(cp), cpchr(cp));
	break;
case REP:
	fprintf(dat, "ylabel('Re %c_%c', 'FontSize', lfs, 'Fontname', fnm);\n",
                fldchr(cp), cpchr(cp));
	break;
case IMP:
	fprintf(dat, "ylabel('Im %c_%c', 'FontSize', lfs, 'Fontname', fnm);\n",
                fldchr(cp), cpchr(cp));
	break;
}
fprintf(dat, "\n");
	if(nbd > 0)
	{
fprintf(dat, "%% Background: shading according to waveguide geometry\n");
fprintf(dat, "maxn = max(n);\n");
fprintf(dat, "minn = min(n);\n");
fprintf(dat, "for i=1:nbd+1\n");
fprintf(dat, "   if i==1 b0=x0;\n");
fprintf(dat, "   else    b0=bd(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nbd+1 b1=x1;\n");
fprintf(dat, "   else        b1=bd(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   pc = 1.0-(ri(i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "   patch([b0 b0 b1 b1], [y1 y0 y0 y1], [pc pc pc], ...\n");
fprintf(dat, "           'LineStyle', 'none');\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	}
fprintf(dat, "%% Add field values\n");
fprintf(dat, "lw = 2.0;\n");
	if(Mlop_Colour==YES) fprintf(dat, "ls = 'b-';\n");
	else                 fprintf(dat, "ls = 'k-';\n");
	for(i=0; i<=nsec-1; ++i)
	{
fprintf(dat, "plot(pos%c%c, %c%c%c%c, ls, 'LineWidth', lw);\n", 
	     dig10(i), dig1(i), 
             fldchr(cp), cpchr(cp), dig10(i), dig1(i));
	}
fprintf(dat, "\n");
fprintf(dat, "%% Include further field section data\n");
fprintf(dat, "%% _______V\n");
fprintf(dat, "\n");
	if(nbd > 0)
	{
fprintf(dat, "%% Indicate permittivity boundaries\n");
fprintf(dat, "for i=1:nbd\n");
fprintf(dat, "   plot([bd(i), bd(i)], [y0 y1], 'k-', 'LineWidth', 0.5);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([x0 x0 x1 x1 x0], ...\n");
fprintf(dat, "     [y1 y0 y0 y1 y1], 'k-', 'LineWidth', alw);\n");
fprintf(dat, "\n");
	}
	return;
}

void mlout_Vsecplot(FILE *dat, Fcomp cp, int nsec, char ext0, char ext1)
{
	int i;
fprintf(dat, "%% Add field values\n");
fprintf(dat, "%% lw = 2.0;\n");
	if(Mlop_Colour==YES) fprintf(dat, "%% ls = 'r-';\n");
	else                 fprintf(dat, "%% ls = 'k:';\n");
	for(i=0; i<=nsec-1; ++i)
	{
fprintf(dat, "plot(pos%c%c%c%c, %c%c%c%c%c%c, ls, 'LineWidth', lw);\n", 
	     dig10(i), dig1(i), ext0, ext1,
             fldchr(cp), cpchr(cp), dig10(i), dig1(i), ext0, ext1);
	}
fprintf(dat, "\n");
	return;
}

/* 1D refractive index profile, plot commands */
void mlout_refindplot(const char *name, FILE *dat, Interval i)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "axis off;\n");
fprintf(dat, "alw = 1.2;\n");
fprintf(dat, "lw = 1.6;\n");
fprintf(dat, "fs = 16;\n");
fprintf(dat, "lfs = 16;\n");
fprintf(dat, "fnm = 'Helvetica';\n");
fprintf(dat, "\n");

fprintf(dat, "x0 = %g;\n", i.x0);
fprintf(dat, "x1 = %g;\n", i.x1);
fprintf(dat, "y0 = minf-(maxf-minf)/20.0;\n");
fprintf(dat, "y1 = maxf+(maxf-minf)/20.0;\n");
fprintf(dat, "%% y0 = minf;\n");
fprintf(dat, "%% y1 = maxf;\n");
fprintf(dat, "\n");

fprintf(dat, "h = axes('Position', [0.15 0.15 0.80 0.45], ...\n");
fprintf(dat, "       'XLim', [x0 x1], ...\n");
fprintf(dat, "       'YLim', [y0 y1], ...\n");
fprintf(dat, "       'FontSize', fs, ...\n");
fprintf(dat, "       'FontName', fnm, ...\n");
fprintf(dat, "       'LineWidth', alw, ...\n");
fprintf(dat, "       'TickDir', 'out');\n");

fprintf(dat, "axes(h);\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");

fprintf(dat, "xlabel(['x [' char(956) 'm]'], 'FontSize', lfs, 'FontName', fnm);\n");
fprintf(dat, "ylabel('n', 'FontSize', lfs, 'Fontname', fnm);\n");
fprintf(dat, "\n");

fprintf(dat, "%% plot waveguide geometry\n");
fprintf(dat, "lw = 2.0;\n");
	if(Mlop_Colour==YES) fprintf(dat, "ls = 'b-';\n");
	else                 fprintf(dat, "ls = 'k-';\n");
fprintf(dat, "maxn = max(n);\n");
fprintf(dat, "minn = min(n);\n");
fprintf(dat, "for i=1:nlX+2\n");
fprintf(dat, "   if i==1 b0=-100.0;\n");
fprintf(dat, "   else    b0=bdX(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nlX+2 b1=100.0;\n");
fprintf(dat, "   else        b1=bdX(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   pc = 1.0-(n(i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "   patch([b0 b0 b1 b1], [y1 y0 y0 y1], [pc pc pc], ...\n");
fprintf(dat, "           'LineStyle', 'none');\n");
fprintf(dat, "   plot([b0 b1], [n(i) n(i)], ls, 'LineWidth', lw);\n"); 
fprintf(dat, "end\n");
fprintf(dat, "\n");

fprintf(dat, "%% Indicate permittivity boundaries\n");
fprintf(dat, "for i=1:nlX+1\n");
fprintf(dat, "   plot([bdX(i), bdX(i)], [y0 y1], 'k-', 'LineWidth', 0.5);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([x0 x0 x1 x1 x0], ...\n");
fprintf(dat, "     [y1 y0 y0 y1 y1], 'k-', 'LineWidth', alw);\n");
fprintf(dat, "\n");

	return;
}

/* 1D complex refractive index profile, plot commands */
void mlout_crefindplot(const char *name, FILE *dat, Interval i)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "axis off;\n");
fprintf(dat, "alw = 1.2;\n");
fprintf(dat, "lw = 1.6;\n");
fprintf(dat, "fs = 16;\n");
fprintf(dat, "lfs = 16;\n");
fprintf(dat, "fnm = 'Helvetica';\n");
fprintf(dat, "\n");

fprintf(dat, "x0 = %g;\n", i.x0);
fprintf(dat, "x1 = %g;\n", i.x1);
fprintf(dat, "y0 = minf-(maxf-minf)/20.0;\n");
fprintf(dat, "y1 = maxf+(maxf-minf)/20.0;\n");
fprintf(dat, "%% y0 = minf;\n");
fprintf(dat, "%% y1 = maxf;\n");
fprintf(dat, "\n");

fprintf(dat, "h = axes('Position', [0.15 0.15 0.80 0.45], ...\n");
fprintf(dat, "       'XLim', [x0 x1], ...\n");
fprintf(dat, "       'YLim', [y0 y1], ...\n");
fprintf(dat, "       'FontSize', fs, ...\n");
fprintf(dat, "       'FontName', fnm, ...\n");
fprintf(dat, "       'LineWidth', alw, ...\n");
fprintf(dat, "       'TickDir', 'out');\n");

fprintf(dat, "axes(h);\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");

fprintf(dat, "xlabel(['x [' char(956) 'm]'], 'FontSize', lfs, 'FontName', fnm);\n");
fprintf(dat, "ylabel('Re n, Im n', 'FontSize', lfs, 'Fontname', fnm);\n");
fprintf(dat, "\n");

fprintf(dat, "%% plot waveguide geometry\n");
fprintf(dat, "lw = 2.0;\n");
	if(Mlop_Colour==YES) { fprintf(dat, "rels = 'b-';\n"); fprintf(dat, "imls = 'k--';\n"); }
	else                 { fprintf(dat, "rels = 'k-';\n"); fprintf(dat, "imls = 'k--';\n"); }
fprintf(dat, "n = abs(n_re + i*n_im);");
fprintf(dat, "maxn = max(n);\n");
fprintf(dat, "minn = min(n);\n");
fprintf(dat, "for i=1:nlX+2\n");
fprintf(dat, "   if i==1 b0=-100.0;\n");
fprintf(dat, "   else    b0=bdX(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nlX+2 b1=100.0;\n");
fprintf(dat, "   else        b1=bdX(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   pc = 1.0-(n(i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "   patch([b0 b0 b1 b1], [y1 y0 y0 y1], [pc pc pc], ...\n");
fprintf(dat, "           'LineStyle', 'none');\n");
fprintf(dat, "   plot([b0 b1], [n_re(i) n_re(i)], rels, 'LineWidth', lw);\n"); 
fprintf(dat, "   plot([b0 b1], [n_im(i) n_im(i)], imls, 'LineWidth', lw);\n"); 
fprintf(dat, "end\n");
fprintf(dat, "\n");

fprintf(dat, "%% Indicate permittivity boundaries\n");
fprintf(dat, "for i=1:nlX+1\n");
fprintf(dat, "   plot([bdX(i), bdX(i)], [y0 y1], 'k-', 'LineWidth', 0.5);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([x0 x0 x1 x1 x0], ...\n");
fprintf(dat, "     [y1 y0 y0 y1 y1], 'k-', 'LineWidth', alw);\n");
fprintf(dat, "\n");

	return;
}

/* 2D refractive index profile, plot commands */
void mlout_refindoverview(const char *name, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  n(x, z)', 'Color', 'w');\n", name);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "lfs = 12;\n");
fprintf(dat, "alw = 1.5;\n");
fprintf(dat, "box on; hold on;\n");
fprintf(dat, "axis([zbeg, zend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(zend-zbeg)/(xend-xbeg) 0.95 0.95]);\n");
fprintf(dat, "xlabel(['z [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Draw refractive index patches\n");
fprintf(dat, "maxn = max(npn);\n");
fprintf(dat, "minn = min(npn);\n");
fprintf(dat, "for i=1:numnp\n");
fprintf(dat, "   pc = 1.0-(npn(i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "   patch([npz0(i) npz0(i) npz1(i) npz1(i)], ...\n");
fprintf(dat, "         [npx1(i) npx0(i) npx0(i) npx1(i)], ...\n");
fprintf(dat, "         [pc pc pc], 'LineStyle', 'none');\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Permittivity boundaries\n");
fprintf(dat, "lw = 1.5; ls = 'k-';\n");
	mlout_geooutline(dat);
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([zbeg zbeg zend zend zbeg], ...\n");
fprintf(dat, "     [xend xbeg xbeg xend xend], 'k-', 'LineWidth', 1.5);\n");
fprintf(dat, "\n");
	return;
}

/* print commands */
void mlout_print(FILE *dat, const char *name, char dev)
{
	if(Mlop_Print != YES) 
	{
		switch(dev)
		{
		case 'p':
fprintf(dat, "%%%% Print figure to file\n");
fprintf(dat, "%% print -dpng -r200 %s \n", name);
			break;
		default:
fprintf(dat, "%%%% Write encapsulated postscript file\n");
fprintf(dat, "%% print -depsc2 -r600 %s\n", name);
fprintf(dat, "\n");
			break;
		}
		return;	
	}
	switch(dev)
	{
	case 'p':
fprintf(dat, "%% Print figure to file\n");
fprintf(dat, "print -dpng -r200 %s \n", name);
		break;
	default:
fprintf(dat, "%% Write encapsulated postscript file\n");
fprintf(dat, "print -depsc2 -r600 %s\n", name);
fprintf(dat, "\n");
		break;
	}
	return;	
}

/* .m file title comment */
void mlout_title(FILE *dat, const char *name, const char *intro)
{
fprintf(dat, "%% %s\n", name);
fprintf(dat, "%% %s\n", intro);
fprintf(dat, "\n");
return;
}

/* matrix of field values */
void mlout_fld(FILE *dat, int npx, int npz, Fcomp cp, Dmatrix fld)
{
        int xi, zi;
        fprintf(dat, "%% field values %c%c\n", fldchr(cp), cpchr(cp));
        fprintf(dat, "%c%c = [\n", fldchr(cp), cpchr(cp));
        for(xi=0; xi<=npx-1; ++xi)
        {
                for(zi=0; zi<=npz-2; ++zi) fprintf(dat, "%g ", fld(xi, zi));
                fprintf(dat, "%g ;\n", fld(xi, npz-1));
        }
        fprintf(dat, "];\n");
        fprintf(dat, "\n");
        return;
}

/* interference plot, general structure */
void mlout_genimage(Fcomp cp, Afo foa, const char *name, const char *desc, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  %s(x, z)', 'Color', 'w');\n", name, desc);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "%% adjust intensity range to actual field\n");
fprintf(dat, "maxf = max(max(%c%c));\n", fldchr(cp), cpchr(cp));
fprintf(dat, "%% minf = min(min(%c%c));\n", fldchr(cp), cpchr(cp));
	switch(foa)
	{
	case MOD: case SQR:
fprintf(dat, "minf = 0.0;\n");
		break;
	case ORG: case REP: case IMP:
fprintf(dat, "minf = -maxf;\n");
		break;
	}
fprintf(dat, "if minf >= maxf\n");
fprintf(dat, "   minf = -1.0;\n");
fprintf(dat, "   maxf =  1.0;\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Intensity image\n");
fprintf(dat, "imagesc(z, x, %c%c, [minf, maxf]);\n", fldchr(cp), cpchr(cp)); 
fprintf(dat, "\n");
	mlout_contours(cp, dat);
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "lfs = 12;\n");
fprintf(dat, "alw = 1.5;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "axis([zbeg, zend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(zend-zbeg)/(xend-xbeg) 0.95 0.95]);\n");
fprintf(dat, "xlabel(['z [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "\n");
	mlout_sketchinit(dat);
	mlout_geooutline(dat);
	return;
}


/* interference plot, general structure, generic */
void mlout_genimage(char xchr, char zchr, Fcomp cp, Afo foa, const char *name, const char *desc, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  %s(%c, %c)', 'Color', 'w');\n", name, desc, xchr, zchr);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "%% adjust intensity range to actual field\n");
fprintf(dat, "maxf = max(max(%c%c));\n", fldchr(cp), cpchr(cp));
fprintf(dat, "%% minf = min(min(%c%c));\n", fldchr(cp), cpchr(cp));
	switch(foa)
	{
	case MOD: case SQR:
fprintf(dat, "minf = 0.0;\n");
		break;
	case ORG: case REP: case IMP:
fprintf(dat, "minf = -maxf;\n");
		break;
	}
fprintf(dat, "if minf >= maxf\n");
fprintf(dat, "   minf = -1.0;\n");
fprintf(dat, "   maxf =  1.0;\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Intensity image\n");
fprintf(dat, "imagesc(%c, %c, %c%c, [minf, maxf]);\n", zchr, xchr, fldchr(cp), cpchr(cp)); 
fprintf(dat, "\n");
	mlout_contours(xchr, zchr, cp, dat);
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "lfs = 12;\n");
fprintf(dat, "alw = 1.5;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "axis([%cbeg, %cend, %cbeg, %cend]);\n", zchr, zchr, xchr, xchr);
fprintf(dat, "set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(%cend-%cbeg)/(%cend-%cbeg) 0.95 0.95]);\n", zchr, zchr, xchr, xchr);
fprintf(dat, "xlabel(['%c [' char(956) 'm]'], 'FontSize', lfs);\n", zchr);
fprintf(dat, "ylabel(['%c [' char(956) 'm]'], 'FontSize', lfs);\n", xchr);
fprintf(dat, "\n");
	mlout_sketchinit(dat);
	mlout_geooutline(dat);
	return;
}

/* interference animation, general structure */
void mlout_genmovie(Fcomp cp, const char *name, const char *desc, FILE *dat, 
                    int ntfr, double dt, double om, double famp)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  %s(x, z, t)', 'Color', 'w');\n", name, desc);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
	mlout_magcolormap(dat);
fprintf(dat, "%% intensity range\n");
fprintf(dat, "maxf = %g;\n", famp);
fprintf(dat, "\n");

fprintf(dat, "%% Loop over one time period\n");
fprintf(dat, "for ti=0:%d\n", ntfr-1);
fprintf(dat, "   t = ti*%.10g;\n", dt);
fprintf(dat, "   fld = %c%c_re*cos(%.10g*t) - %c%c_im*sin(%.10g*t);\n",
             fldchr(cp), cpchr(cp), om, fldchr(cp), cpchr(cp), om); 

fprintf(dat, "   %% Intensity image\n");
fprintf(dat, "   imagesc(z, x, fld, [-maxf, maxf]);\n");
fprintf(dat, "\n");
	mlout_contours(dat);
fprintf(dat, "   %% Adjust axes and labels\n");
fprintf(dat, "   lfs = 12;\n");
fprintf(dat, "   alw = 1.5;\n");
fprintf(dat, "   box on;\n");
fprintf(dat, "   axis([zbeg, zend, xbeg, xend]);\n");
fprintf(dat, "   set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(zend-zbeg)/(xend-xbeg) 0.95 0.95]);\n");
fprintf(dat, "   xlabel(['z [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "   ylabel(['x [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "\n");
	mlout_sketchinit(dat);
	mlout_geooutline(dat);
fprintf(dat, "   %% show time-position\n");
fprintf(dat, "   text(zend-(zend-zbeg)/30, xbeg+(xend-xbeg)/30, ...\n");
fprintf(dat, "      sprintf('t = %%#4.2f fs', t), 'Color', 'w', ...\n");
fprintf(dat, "      'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "      'VerticalAlignment', 'bottom', 'FontSize', 12);\n");  
fprintf(dat, "   \n");
fprintf(dat, "   if ti == 0\n");
fprintf(dat, "      %% reserve memory\n");
fprintf(dat, "      M = moviein(%d);\n", ntfr);     
fprintf(dat, "   end\n");
fprintf(dat, "   %% store frames\n");
fprintf(dat, "   M(:,ti+1) = getframe;\n");  
fprintf(dat, "   hold off;\n");
	if(Mlop_Print != YES) 
	{
fprintf(dat, "   %%%% Print frame to file\n");
fprintf(dat, "   %% print('-dpng', '-r75', sprintf('%s_f%%d',ti));\n", name);
	}
	else
	{
fprintf(dat, "   %% Print figure to file\n");
fprintf(dat, "   print('-dpng', '-r75', sprintf('%s_f%%d',ti));\n", name);
	}
fprintf(dat, "\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	if(Mlop_Print != YES) 
	{
fprintf(dat, "%%%% Save movie as a video file\n");
fprintf(dat, "%%movie2avi(M, '%s', ...\n", name);
fprintf(dat, "%%   'fps', 15, ...\n");
fprintf(dat, "%%   'compression', 'none', ...\n");
fprintf(dat, "%%   'quality', 100, ...\n");
fprintf(dat, "%%   'videoname', '%s(x, z, t)');\n", desc);
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "movie(M, 20)\n");         
	}
	else
	{
fprintf(dat, "%% Save movie as a video file\n");
fprintf(dat, "movie2avi(M, '%s', ...\n", name);
fprintf(dat, "   'fps', 15, ...\n");
fprintf(dat, "   'compression', 'none', ...\n");
fprintf(dat, "   'quality', 100, ...\n");
fprintf(dat, "   'videoname', '%s(x, z, t)');\n", desc);
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "%% movie(M, 20)\n");         
	}
fprintf(dat, "\n");
	return;
}

/* interference animation, general structure, generic coordinates */
void mlout_genmovie(char xchr, char zchr, Fcomp cp, 
                    const char *name, const char *desc, FILE *dat, 
                    int ntfr, double dt, double om, double famp)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  %s(%c, %c, t)', 'Color', 'w');\n", name, desc, zchr, xchr);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
	mlout_magcolormap(dat);
fprintf(dat, "%% intensity range\n");
fprintf(dat, "maxf = %g;\n", famp);
fprintf(dat, "\n");

fprintf(dat, "%% Loop over one time period\n");
fprintf(dat, "for ti=0:%d\n", ntfr-1);
fprintf(dat, "   t = ti*%.10g;\n", dt);
fprintf(dat, "   fld = %c%c_re*cos(%.10g*t) - %c%c_im*sin(%.10g*t);\n",
             fldchr(cp), cpchr(cp), om, fldchr(cp), cpchr(cp), om); 

fprintf(dat, "   %% Intensity image\n");
fprintf(dat, "   imagesc(%c, %c, fld, [-maxf, maxf]);\n", zchr, xchr);
fprintf(dat, "\n");
	mlout_contours(xchr, zchr, dat);
fprintf(dat, "   %% Adjust axes and labels\n");
fprintf(dat, "   lfs = 12;\n");
fprintf(dat, "   alw = 1.5;\n");
fprintf(dat, "   box on;\n");
fprintf(dat, "   axis([%cbeg, %cend, %cbeg, %cend]);\n", zchr, zchr, xchr, xchr);
fprintf(dat, "   set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(%cend-%cbeg)/(%cend-%cbeg) 0.95 0.95]);\n", zchr, zchr, xchr, xchr);
fprintf(dat, "   xlabel(['%c [' char(956) 'm]'], 'FontSize', lfs);\n", zchr);
fprintf(dat, "   ylabel(['%c [' char(956) 'm]'], 'FontSize', lfs);\n", xchr);
fprintf(dat, "\n");
	mlout_sketchinit(dat);
	mlout_geooutline(dat);
fprintf(dat, "   %% show time-position\n");
fprintf(dat, "   text(%cend-(%cend-%cbeg)/30, %cbeg+(%cend-%cbeg)/30, ...\n", zchr, zchr, zchr, xchr, xchr, xchr);
fprintf(dat, "      sprintf('t = %%#4.2f fs', t), 'Color', 'w', ...\n");
fprintf(dat, "      'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "      'VerticalAlignment', 'bottom', 'FontSize', 12);\n");  
fprintf(dat, "   \n");
fprintf(dat, "   if ti == 0\n");
fprintf(dat, "      %% reserve memory\n");
fprintf(dat, "      M = moviein(%d);\n", ntfr);     
fprintf(dat, "   end\n");
fprintf(dat, "   %% store frames\n");
fprintf(dat, "   M(:,ti+1) = getframe;\n");  
fprintf(dat, "   hold off;\n");
	if(Mlop_Print != YES) 
	{
fprintf(dat, "   %%%% Print frame to file\n");
fprintf(dat, "   %% print('-dpng', '-r75', sprintf('%s_f%%d',ti));\n", name);
	}
	else
	{
fprintf(dat, "   %% Print figure to file\n");
fprintf(dat, "   print('-dpng', '-r75', sprintf('%s_f%%d',ti));\n", name);
	}
fprintf(dat, "\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	if(Mlop_Print != YES) 
	{
fprintf(dat, "%%%% Save movie as a video file\n");
fprintf(dat, "%%movie2avi(M, '%s', ...\n", name);
fprintf(dat, "%%   'fps', 15, ...\n");
fprintf(dat, "%%   'compression', 'none', ...\n");
fprintf(dat, "%%   'quality', 100, ...\n");
fprintf(dat, "%%   'videoname', '%s(x, z, t)');\n", desc);
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "movie(M, 20)\n");         
	}
	else
	{
fprintf(dat, "%% Save movie as a video file\n");
fprintf(dat, "movie2avi(M, '%s', ...\n", name);
fprintf(dat, "   'fps', 15, ...\n");
fprintf(dat, "   'compression', 'none', ...\n");
fprintf(dat, "   'quality', 100, ...\n");
fprintf(dat, "   'videoname', '%s(x, z, t)');\n", desc);
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "%% movie(M, 20)\n");         
	}
fprintf(dat, "\n");
	return;
}

/* fancy interference plot */
void mlout_fancy(const char *name, const char *desc, FILE *dat, 
                 Fcomp cp, int nc)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  %s(x, z)', 'Color', 'w');\n", name, desc);
fprintf(dat, "\n");
fprintf(dat, "%% Adjust field aspect ratio to cross section extension\n");
fprintf(dat, "s = 1; %% set to -1 to invert the plot \n");
fprintf(dat, "fld = %c%c;\n", fldchr(cp), cpchr(cp));
fprintf(dat, "xyzfac = min([zend-zbeg, xend-xbeg])*0.3;\n");
fprintf(dat, "maxf = max(max(fld));\n");
fprintf(dat, "minf = min(min(fld));\n");
fprintf(dat, "df   = maxf-minf;\n");
fprintf(dat, "fld  =  s*fld/df*xyzfac;\n");
fprintf(dat, "minf = s*minf/df*xyzfac;\n");
fprintf(dat, "maxf = s*maxf/df*xyzfac;\n");
fprintf(dat, "if minf > maxf tf=minf; minf=maxf; maxf=tf; end\n");
fprintf(dat, "\n");
	mlout_magcolormap(dat);
fprintf(dat, "\n");
fprintf(dat, "%% Surface plot \n");
fprintf(dat, "sf = surf(z, x, fld);\n");
fprintf(dat, "set(sf, 'LineStyle', 'none');\n");
fprintf(dat, "set(sf, 'FaceColor', 'interp');\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "axis([zbeg, zend, xbeg, xend, minf, maxf]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1], ...\n");
fprintf(dat, "         'Projection', 'Perspective', ...\n");
fprintf(dat, "         'Units', 'normalized', ...\n");
fprintf(dat, "         'Position', [0, 0, 1, 1]);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "grid off;\n");
fprintf(dat, "view(30,25);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Lighting\n");
fprintf(dat, "lighting phong;\n");
fprintf(dat, "light('Position', [(zend+zbeg)/2, xbeg-2*(xend-xbeg), maxf+2*(maxf-minf)]);\n");
fprintf(dat, "light('Position', [(zend+zbeg)/2, xend+2*(xend-xbeg), maxf+2*(maxf-minf)]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries \n");
fprintf(dat, "hold on;\n");
fprintf(dat, "lw = 3.0;\n");
fprintf(dat, "lc = 'k';\n");
fprintf(dat, "lift = s*df/200.0;\n");
fprintf(dat, "for j=1:nc\n");
fprintf(dat, "   zp = linspace(bdz0(j), bdz1(j), bdnp(j));\n");
fprintf(dat, "   xp = linspace(bdx0(j), bdx1(j), bdnp(j));\n");
fprintf(dat, "   f = bdf%c%c(bdsi(j):bdsi(j)+bdnp(j)-1);\n", 
             fldchr(cp), cpchr(cp));
fprintf(dat, "   line(zp, xp, s*(f+lift)/df*xyzfac, 'LineWidth', lw, 'Color', lc);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* interference animation, fancy style */
void mlout_genfmovie(Fcomp cp, const char *name, const char *desc, FILE *dat,
                     int nc, int ntfr, double dt, double om, double famp)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  %s(x, z, t)', 'Color', 'w');\n", name, desc);
fprintf(dat, "\n");
fprintf(dat, "%% Adjust field aspect ratio to cross section extension\n");
fprintf(dat, "s = 1; %% set to -1 to invert the plot \n");
fprintf(dat, "fld_re = %c%c_re;\n", fldchr(cp), cpchr(cp));
fprintf(dat, "fld_im = %c%c_im;\n", fldchr(cp), cpchr(cp));
fprintf(dat, "xyzfac = min([zend-zbeg, xend-xbeg])*0.3;\n");
fprintf(dat, "maxf = %g;\n", famp);
fprintf(dat, "minf = %g;\n", -famp);
fprintf(dat, "df   = maxf-minf;\n");
fprintf(dat, "fld_re = s*fld_re/df*xyzfac;\n");
fprintf(dat, "fld_im = s*fld_im/df*xyzfac;\n");
fprintf(dat, "minf = s*minf/df*xyzfac;\n");
fprintf(dat, "maxf = s*maxf/df*xyzfac;\n");
fprintf(dat, "if minf > maxf tf=minf; minf=maxf; maxf=tf; end\n");
fprintf(dat, "\n");
	mlout_magcolormap(dat);
fprintf(dat, "%% Loop over one time period\n");
fprintf(dat, "for ti=0:%d\n", ntfr-1);
fprintf(dat, "   t = ti*%.10g;\n", dt);
fprintf(dat, "   fld = fld_re*cos(%.10g*t) - fld_im*sin(%.10g*t);\n", om, om); 
fprintf(dat, "\n");
fprintf(dat, "%% Surface plot \n");
fprintf(dat, "sf = surf(z, x, fld);\n");
fprintf(dat, "set(sf, 'LineStyle', 'none');\n");
fprintf(dat, "set(sf, 'FaceColor', 'interp');\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "axis([zbeg, zend, xbeg, xend, minf, maxf]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1], ...\n");
fprintf(dat, "         'Projection', 'Perspective', ...\n");
fprintf(dat, "         'Units', 'normalized', ...\n");
fprintf(dat, "         'Position', [0, 0, 1, 1]);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "grid off;\n");
fprintf(dat, "view(30,25);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Lighting\n");
fprintf(dat, "lighting phong;\n");
fprintf(dat, "light('Position', [(zend+zbeg)/2, xbeg-2*(xend-xbeg), maxf+2*(maxf-minf)]);\n");
fprintf(dat, "light('Position', [(zend+zbeg)/2, xend+2*(xend-xbeg), maxf+2*(maxf-minf)]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries \n");
fprintf(dat, "hold on;\n");
fprintf(dat, "lw = 3.0;\n");
fprintf(dat, "lc = 'k';\n");
fprintf(dat, "lift = s*df/200.0;\n");
fprintf(dat, "for j=1:nc\n");
fprintf(dat, "   zp = linspace(bdz0(j), bdz1(j), bdnp(j));\n");
fprintf(dat, "   xp = linspace(bdx0(j), bdx1(j), bdnp(j));\n");
fprintf(dat, "   f_re = bdf%c%c_re(bdsi(j):bdsi(j)+bdnp(j)-1);\n",
             fldchr(cp), cpchr(cp));
fprintf(dat, "   f_im = bdf%c%c_im(bdsi(j):bdsi(j)+bdnp(j)-1);\n",
             fldchr(cp), cpchr(cp));
fprintf(dat, "   sfld = f_re*cos(%.10g*t) - f_im*sin(%.10g*t);\n", om, om); 
fprintf(dat, "   line(zp, xp, s*(sfld+lift)/df*xyzfac, 'LineWidth', lw, 'Color', lc);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "   if ti == 0\n");
fprintf(dat, "      %% reserve memory\n");
fprintf(dat, "      M = moviein(%d);\n", ntfr);     
fprintf(dat, "   end\n");
fprintf(dat, "   %% store frames\n");
fprintf(dat, "   M(:,ti+1) = getframe;\n");  
fprintf(dat, "   hold off;\n");
	if(Mlop_Print != YES) 
	{
fprintf(dat, "   %%%% Print frame to file\n");
fprintf(dat, "   %% print('-dpng', '-r75', sprintf('%s_f%%d',ti));\n", name);
	}
	else
	{
fprintf(dat, "   %% Print figure to file\n");
fprintf(dat, "   print('-dpng', '-r75', sprintf('%s_f%%d',ti));\n", name);
	}
fprintf(dat, "\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	if(Mlop_Print != YES) 
	{
fprintf(dat, "%%%% Save movie as a video file\n");
fprintf(dat, "%%movie2avi(M, '%s', ...\n", name);
fprintf(dat, "%%   'fps', 15, ...\n");
fprintf(dat, "%%   'compression', 'none', ...\n");
fprintf(dat, "%%   'quality', 100, ...\n");
fprintf(dat, "%%   'videoname', '%s(x, z, t)');\n", desc);
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "movie(M, 20)\n");         
	}
	else
	{
fprintf(dat, "%% Save movie as a video file\n");
fprintf(dat, "movie2avi(M, '%s', ...\n", name);
fprintf(dat, "   'fps', 15, ...\n");
fprintf(dat, "   'compression', 'none', ...\n");
fprintf(dat, "   'quality', 100, ...\n");
fprintf(dat, "   'videoname', '%s(x, z, t)');\n", desc);
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "%% movie(M, 20)\n");         
	}
fprintf(dat, "\n");
	return;
}

/* ------------------------------------------------------------------------ */               
/* long three segment coupler: geometry for interference visualization */
void mlout_l3scgeoxz_2(FILE *dat, Waveguide wg1, Waveguide wg2, Waveguide wg3,
                       double l, double xbeg, double xend) 
{
	int i, j;
	double r;
	fprintf(dat, "%% Waveguide geometry\n");
	fprintf(dat, "nlX1 = %d;\n", wg1.nx);
	fprintf(dat, "bdX1 = [\n");
	for(j=0; j<=wg1.nx; ++j) fprintf(dat, "%g ...\n", wg1.hx(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Refractive indices\n");
	fprintf(dat, "n1 = [\n");
	for(j=0; j<=wg1.nx+1; ++j) fprintf(dat, "%g ...\n", wg1.n(j));
	fprintf(dat, "];\n");
	fprintf(dat, "nlX2 = %d;\n", wg2.nx);
	fprintf(dat, "bdX2 = [\n");
	for(j=0; j<=wg2.nx; ++j) fprintf(dat, "%g ...\n", wg2.hx(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Refractive indices\n");
	fprintf(dat, "n2 = [\n");
	for(j=0; j<=wg2.nx+1; ++j) fprintf(dat, "%g ...\n", wg2.n(j));
	fprintf(dat, "];\n");
	fprintf(dat, "nlX3 = %d;\n", wg3.nx);
	fprintf(dat, "bdX3 = [\n");
	for(j=0; j<=wg3.nx; ++j) fprintf(dat, "%g ...\n", wg3.hx(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Refractive indices\n");
	fprintf(dat, "n3 = [\n");
	for(j=0; j<=wg3.nx+1; ++j) fprintf(dat, "%g ...\n", wg3.n(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	fprintf(dat, "%% Vertical discontinuities at z=0\n");
	Dvector xpos;
	Dvector beg;
	Dvector end;
	int x1, x2, p;
	double x1p, x2p;
	x1 = 0; x2 = 0; xpos.strip();
	while(x1 <= wg1.nx || x2 <= wg2.nx)
	{
		if(x1<=wg1.nx) x1p = wg1.hx(x1); else x1p=1.0e+300;
		if(x2<=wg2.nx) x2p = wg2.hx(x2); else x2p=1.0e+300;
		if(x1p < x2p)
		{
			xpos.append(x1p);
			++x1;
		}
		else
		{
			xpos.append(x2p);
			++x2;
		}
	}
	p = xpos.nel;
	for(i=-1; i<=p-1; ++i)
	{
		if(i==-1) x1p = xbeg; 
		else      x1p = xpos(i);
		if(i==p-1) x2p = xend;
		else       x2p = xpos(i+1);
		if(fabs(x1p-x2p) >= 1.0e-10)
		{
			r = 0.5*(x1p+x2p);
			if(fabs(wg1.eps(r)-wg2.eps(r))>1.0e-10)
			{
				beg.append(x1p);
				end.append(x2p);
			} 
		}
	}
	j = beg.nel;
	fprintf(dat, "nvd0 = %d;\n", j);
	fprintf(dat, "vd0beg = [\n");
	for(i=0; i<=j-1; ++i) fprintf(dat, "%g ...\n", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vd0end = [\n");
	for(i=0; i<=j-1; ++i) fprintf(dat, "%g ...\n", end(i));
	fprintf(dat, "];\n");
	beg.strip();
	end.strip();

	x1 = 0; x2 = 0; xpos.strip(); 
	while(x1 <= wg2.nx || x2 <= wg3.nx)
	{
		if(x1<=wg2.nx) x1p = wg2.hx(x1); else x1p=1.0e+300;
		if(x2<=wg3.nx) x2p = wg3.hx(x2); else x2p=1.0e+300;
		if(x1p < x2p)
		{
			xpos.append(x1p);
			++x1;
		}
		else
		{
			xpos.append(x2p);
			++x2;
		}
	}
	p = xpos.nel;
	for(i=-1; i<=p-1; ++i)
	{
		if(i==-1) x1p = xbeg; 
		else      x1p = xpos(i);
		if(i==p-1) x2p = xend;
		else       x2p = xpos(i+1);
		if(fabs(x1p-x2p) >= 1.0e-10)
		{
			r = 0.5*(x1p+x2p);
			if(fabs(wg2.eps(r)-wg3.eps(r))>1.0e-10)
			{
				beg.append(x1p);
				end.append(x2p);
			} 
			
		}
	}
	j = beg.nel;
	fprintf(dat, "nvdl = %d;\n", j);
	fprintf(dat, "vdlbeg = [\n");
	for(i=0; i<=j-1; ++i) fprintf(dat, "%g ...\n", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdlend = [\n");
	for(i=0; i<=j-1; ++i) fprintf(dat, "%g ...\n", end(i));
	fprintf(dat, "];\n");

	fprintf(dat, "\n");
	fprintf(dat, "%% Length of the coupling segment\n");
	fprintf(dat, "L = %g;\n", l);
	fprintf(dat, "\n");
	return;
}

/* ------------------------------------------------------------------------ */               
/* interference intensity image, long three segment coupler */
void mlout_1Dl3scimage(const char *name, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
	mlout_lincolormap(dat);
fprintf(dat, "\n");
fprintf(dat, "%% adjust intensity range to actual field\n");
fprintf(dat, "maxf = max(max(Sz));\n");
fprintf(dat, "minf = min(min(Sz));\n");
fprintf(dat, "if minf >= maxf\n");
fprintf(dat, "   minf = minf-1.0;\n");
fprintf(dat, "   maxf = minf+2.0;\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Intensity image\n");
fprintf(dat, "imagesc(z, x, Sz, [minf, maxf]);\n"); 
fprintf(dat, "\n");
	mlout_contours(SZ, dat);
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "lfs = 12;\n");
fprintf(dat, "alw = 1.5;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "axis([zbeg, zend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(zend-zbeg)/(xend-xbeg) 0.95 0.95]);\n");
fprintf(dat, "xlabel(['z [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "\n");
	mlout_sketchinit(dat);
fprintf(dat, "for i=1:nlX1+1\n");
fprintf(dat, "   xp=bdX1(i); \n");
fprintf(dat, "   if n1(i)~=n1(i+1)\n");
fprintf(dat, "      plot([zbeg, 0], [xp, xp], ls, 'LineWidth', lw);\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for i=1:nlX2+1\n");
fprintf(dat, "   xp=bdX2(i); \n");
fprintf(dat, "   if n2(i)~=n2(i+1)\n");
fprintf(dat, "      plot([0, L], [xp, xp], ls, 'LineWidth', lw);\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for i=1:nlX3+1\n");
fprintf(dat, "   xp=bdX3(i); \n");
fprintf(dat, "   if n3(i)~=n3(i+1)\n");
fprintf(dat, "      plot([L, zend], [xp, xp], ls, 'LineWidth', lw);\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for i=1:nvd0\n");
fprintf(dat, "   plot([0, 0], [vd0beg(i), vd0end(i)], ls, 'LineWidth', lw);\n");
fprintf(dat, "end\n");
fprintf(dat, "for i=1:nvdl\n");
fprintf(dat, "   plot([L, L], [vdlbeg(i), vdlend(i)], ls, 'LineWidth', lw);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* ------------------------------------------------------------------------ */ 

/* field viewer .m-file: top, title and globals */
void mlout_viewertop(FILE *dat, const char *name, Polarization p, double lambda)
{
fprintf(dat, "%% %s.m\n", name);
fprintf(dat, "%% Metric field solution, viewer GUI\n");
fprintf(dat, "function %s\n", name);
fprintf(dat, "\n");
fprintf(dat, "close all; clear global;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Polarisation\n");
fprintf(dat, "pol = 'T%c';\n", polCHR(p));
fprintf(dat, "%% Wavelength\n");
fprintf(dat, "lambda = %.10g;\n", lambda);
fprintf(dat, "\n");
fprintf(dat, "%% Field and structure data\n");
fprintf(dat, "global nvd vdbeg vdend vdpos nhd hdbeg hdend hdpos n;\n");
fprintf(dat, "global ncd cdrad cdvp cdhp;\n");
fprintf(dat, "global xbeg xend numx x zbeg zend numz z;\n");
fprintf(dat, "global Ex_re Ex_im Ey_re Ey_im Ez_re Ez_im;\n");
fprintf(dat, "global Hx_re Hx_im Hy_re Hy_im Hz_re Hz_im;\n");
fprintf(dat, "\n");
return;
}
void mlout_viewertop(FILE *dat, const char *name, Polarization p, double lambda, int vM)
{
	mlout_viewertop(dat, name, 'x', 'z', p, lambda, vM);
	return;
}
void mlout_viewertop(FILE *dat, const char *name, char xchr, char zchr, Polarization p, double lambda, int vM)
{
fprintf(dat, "%% %s.m\n", name);
fprintf(dat, "%% Metric field solution, viewer GUI\n");
fprintf(dat, "function %s\n", name);
fprintf(dat, "\n");
fprintf(dat, "close all; clear global;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Polarisation\n");
if(vM) 
fprintf(dat, "pol = 'vM';\n");
else 
fprintf(dat, "pol = 'T%c';\n", polCHR(p));
fprintf(dat, "%% Wavelength\n");
fprintf(dat, "lambda = %.10g;\n", lambda);
fprintf(dat, "\n");
fprintf(dat, "%% Field and structure data\n");
fprintf(dat, "global nvd vdbeg vdend vdpos nhd hdbeg hdend hdpos n;\n");
fprintf(dat, "global ncd cdrad cdvp cdhp;\n");
fprintf(dat, "global %cbeg %cend num%c %c %cbeg %cend num%c %c;\n", xchr, xchr, xchr, xchr, zchr, zchr, zchr, zchr);
fprintf(dat, "global Ex_re Ex_im Ey_re Ey_im Ez_re Ez_im;\n");
fprintf(dat, "global Hx_re Hx_im Hy_re Hy_im Hz_re Hz_im;\n");
fprintf(dat, "\n");
return;
}

/* matrix of refractive index values */
void mlout_n(FILE *dat, int npx, int npz, Dmatrix n)
{
        int xi, zi;
        fprintf(dat, "%% Refractive index levels \n");
        fprintf(dat, "n = [\n");
        for(xi=0; xi<=npx-1; ++xi)
        {
                for(zi=0; zi<=npz-2; ++zi) fprintf(dat, "%g ", n(xi, zi));
                fprintf(dat, "%g ;\n", n(xi, npz-1));
        }
        fprintf(dat, "];\n");
        fprintf(dat, "\n");
        return;
}

/* matrix of field values, zero field */
void mlout_0fld(FILE *dat, Fcomp cp)
{
	mlout_0fld(dat, cp, 'x', 'z');
        return;
}
void mlout_0fldxy(FILE *dat, Fcomp cp)
{
	mlout_0fld(dat, cp, 'x', 'y');
        return;
}
void mlout_0fld(FILE *dat, Fcomp cp, char xchr, char zchr)
{
        fprintf(dat, "%% Field values %c%c\n", fldchr(cp), cpchr(cp));
        fprintf(dat, "%c%c = zeros(num%c, num%c);\n", fldchr(cp), cpchr(cp), xchr, zchr);
        fprintf(dat, "\n");
        return;
}

/* field viewer .m-file: viewer commands */
void mlout_fldviewer(FILE *dat, const char *name)
{
	mlout_fldviewer(dat, name, 'x', 'z');
}

/* field viewer .m-file, generic coordinates: viewer commands */
void mlout_fldviewer(FILE *dat, const char *name, char xchr, char zchr)
{
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% Viewer commands \n");
fprintf(dat, "\n");
fprintf(dat, "clear Ex; clear Ey; clear Ez; clear Hx; clear Hy; clear Hz;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Poynting vector\n");
fprintf(dat, "global Sx Sy Sz;\n");
fprintf(dat, "Sx=0.5*real((Ey_re+i*Ey_im).*(Hz_re-i*Hz_im)-(Ez_re+i*Ez_im).*(Hy_re-i*Hy_im)); \n");
fprintf(dat, "Sy=0.5*real((Ez_re+i*Ez_im).*(Hx_re-i*Hx_im)-(Ex_re+i*Ex_im).*(Hz_re-i*Hz_im)); \n");
fprintf(dat, "Sz=0.5*real((Ex_re+i*Ex_im).*(Hy_re-i*Hy_im)-(Ey_re+i*Ey_im).*(Hx_re-i*Hx_im)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Energy density of the optical field\n");
fprintf(dat, "global W;\n");
fprintf(dat, "mu0 = 1.2566370614E+3; ep0 = 8.854187817E-3;\n");
fprintf(dat, "epsr = sign(n).*n.*n;\n");
fprintf(dat, "W = 0.25*(mu0*(Hx_re.*Hx_re+Hx_im.*Hx_im+Hy_re.*Hy_re+Hy_im.*Hy_im ...\n");
fprintf(dat, "              +Hz_re.*Hz_re+Hz_im.*Hz_im ) ...\n");
fprintf(dat, "         +ep0.*epsr.*(Ex_re.*Ex_re+Ex_im.*Ex_im+Ey_re.*Ey_re+Ey_im.*Ey_im ...\n");
fprintf(dat, "                     +Ez_re.*Ez_re+Ez_im.*Ez_im ));\n");
fprintf(dat, "\n");
fprintf(dat, "%% Maximum field levels\n");
fprintf(dat, "global max_E max_H max_S max_W min_n max_n;\n");
fprintf(dat, "max_E = max(max(Ex_re.*Ex_re+Ex_im.*Ex_im));\n");
fprintf(dat, "m = max(max(Ey_re.*Ey_re+Ey_im.*Ey_im));\n");
fprintf(dat, "if m>max_E max_E = m; end\n");
fprintf(dat, "m = max(max(Ez_re.*Ez_re+Ez_im.*Ez_im));\n");
fprintf(dat, "if m>max_E max_E = m; end\n");
fprintf(dat, "max_E = sqrt(max_E);\n");
fprintf(dat, "max_H = max(max(Hx_re.*Hx_re+Hx_im.*Hx_im));\n");
fprintf(dat, "m = max(max(Hy_re.*Hy_re+Hy_im.*Hy_im));\n");
fprintf(dat, "if m>max_H max_H = m; end\n");
fprintf(dat, "m = max(max(Hz_re.*Hz_re+Hz_im.*Hz_im));\n");
fprintf(dat, "if m>max_H max_H = m; end\n");
fprintf(dat, "max_H = sqrt(max_H);\n");
fprintf(dat, "max_S = max(max(abs(Sx)));\n");
fprintf(dat, "m = max(max(abs(Sy)));\n");
fprintf(dat, "if m>max_S max_S = m; end\n");
fprintf(dat, "m = max(max(abs(Sz)));\n");
fprintf(dat, "if m>max_S max_S = m; end\n");
fprintf(dat, "max_W = max(max(abs(W)));\n");
fprintf(dat, "max_n = max(max(n));\n");
fprintf(dat, "min_n = min(min(n));\n");
fprintf(dat, "if min_n >= max_n max_n = min_n+1; end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Linear colormap\n");
fprintf(dat, "global lin_blue;\n");
fprintf(dat, "mi = 32; mlv = 0.2; lin_blue = zeros(64, 3);\n");
fprintf(dat, "lin_blue(1:mi, 3) = linspace(0, 1,   mi)';\n");
fprintf(dat, "lin_blue(1:mi, 2) = linspace(0, mlv, mi)';\n");
fprintf(dat, "lin_blue(1:mi, 1) = linspace(0, mlv, mi)';\n");
fprintf(dat, "lin_blue(mi+1:64, 3) = 1;\n");
fprintf(dat, "lin_blue(mi+1:64, 2) = linspace(mlv, 1, 64-mi)';\n");
fprintf(dat, "lin_blue(mi+1:64, 1) = linspace(mlv, 1, 64-mi)';\n");
fprintf(dat, "\n");
fprintf(dat, "%% Colormap: magnified zero level\n");
fprintf(dat, "global mag_blue;\n");
fprintf(dat, "mag_blue = zeros(256, 3);\n");
fprintf(dat, "cpd = 0.02; cphb= 0.75; cpho= 0.25; cpa = atan(256*cpd/2);\n");
fprintf(dat, "for idx = 0:255\n");
fprintf(dat, "   t = 0.5*(atan(cpd*(idx-256/2))/cpa+1);\n");
fprintf(dat, "   mag_blue(idx+1, 3) = (2-4*cphb)*t*t+(4*cphb-1)*t;\n");
fprintf(dat, "   mag_blue(idx+1, 1) = (2-4*cpho)*t*t+(4*cpho-1)*t;\n");
fprintf(dat, "   mag_blue(idx+1, 2) = mag_blue(idx+1, 1);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% A few style variables ...\n");
fprintf(dat, "global col_fbg col_txt col_abg col_atx col_but col_btx fntsz;\n");
fprintf(dat, "col_fbg = [0.77 0.80 0.87]; col_txt = [0.00 0.00 0.00];\n");
fprintf(dat, "col_abg = [1.00 1.00 1.00]; col_atx = [0.00 0.00 0.00];\n");
fprintf(dat, "col_but = [0.00 0.47 0.48]; col_btx = [1.00 1.00 1.00];\n");
fprintf(dat, "fntsz = 10;\n");
fprintf(dat, "\n");
fprintf(dat, "%% The viewer window\n");
fprintf(dat, "global fgr;\n");
fprintf(dat, "fgr = figure('NumberTitle', 'off', 'Name', '%s', 'Color', col_fbg, ...\n", name);
fprintf(dat, "             'Units', 'pixels', 'Menubar', 'none');\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust figure extension, respect aspect ratio of the data as far as possible,\n");
fprintf(dat, "%% Select figure position\n");
fprintf(dat, "dx = %cend-%cbeg; dz = %cend-%cbeg;\n", xchr, xchr, zchr, zchr);
fprintf(dat, "if dx<0.25*dz dx = 0.25*dz; end\n");
fprintf(dat, "if dz<0.25*dx dz = 0.25*dx; end\n");
fprintf(dat, "axw = round(sqrt(144400*dz/dx)); axh = round(axw*dx/dz);\n");
fprintf(dat, "global csf_xaw csf_zaw;\n");
fprintf(dat, "csf_zaw = axw; csf_xaw = axh;\n");
fprintf(dat, "fgw = axw+80+180; fgh = axh+2*40;\n");
fprintf(dat, "axw = axw/fgw; axh = axh/fgh;\n");
fprintf(dat, "axpx = 60/fgw; axpy = 40/fgh;\n");
fprintf(dat, "cmpx = axpx+axw+7/fgw; cmpw = 13/fgw;\n");
fprintf(dat, "bupx = cmpx+cmpw+45/fgw; buw = 1.0-10/fgw-bupx; busep = 2/fgw;\n");
fprintf(dat, "set(0, 'Units', 'pixels');\n");
fprintf(dat, "scrsz = get(0, 'ScreenSize');\n");
fprintf(dat, "if fgw > scrsz(3)*0.8 \n");
fprintf(dat, "   s = scrsz(3)*0.8/fgw; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "if fgh > scrsz(4)*0.8 \n");
fprintf(dat, "   s = scrsz(4)*0.8/fgh; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "fgpx = round(scrsz(1)+scrsz(3)*0.05); fgpy = round(scrsz(2)+scrsz(4)*0.95)-fgh;\n");
fprintf(dat, "set(fgr, 'Position', [fgpx fgpy fgw fgh]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Plot container: global state \n");
fprintf(dat, "global pl_fld pl_fldre pl_fldim pl_tp pl_min pl_max pl_amin pl_amax;\n");
fprintf(dat, "global pl_cp pl_img pl_afo pl_tih an_nump;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Initial plot: basic field component \n");
fprintf(dat, "if strcmp(pol, 'TE') || strcmp(pol, 'vM') \n");
fprintf(dat, "   pl_tp = 0.0; pl_cp = 'E_y'; pl_afo = 'f';\n");
fprintf(dat, "   pl_fldre = Ey_re; pl_fldim = Ey_im; pl_min = -max_E; pl_max =  max_E;\n");
fprintf(dat, "else\n");
fprintf(dat, "   pl_tp = 0.0; pl_cp = 'H_y'; pl_afo = 'f';\n");
fprintf(dat, "   pl_fldre = Hy_re; pl_fldim = Hy_im; pl_min = -max_H; pl_max =  max_H;\n");
fprintf(dat, "end\n");
fprintf(dat, "an_nump = 10;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Refine the plot window\n");
fprintf(dat, "global axs;\n");
fprintf(dat, "axs = axes;\n");
fprintf(dat, "pl_fld = pl_fldre*cos(2*pi*pl_tp)-pl_fldim*sin(2*pi*pl_tp);\n");
fprintf(dat, "pl_amin = min(min(pl_fld)); pl_amax = max(max(pl_fld));\n");
fprintf(dat, "pl_img = imagesc(%c, %c, pl_fld, [pl_min, pl_max]);\n", zchr, xchr);
fprintf(dat, "set(axs, 'Position', [axpx axpy axw axh], 'Units', 'normalized', ...\n");
fprintf(dat, "         'YDir', 'normal', 'ButtonDownFcn', {@clbf_crossh}, ...\n");
fprintf(dat, "         'XLim', [%cbeg, %cend], 'YLim', [%cbeg, %cend], ...\n", zchr, zchr, xchr, xchr);
fprintf(dat, "         'TickDir', 'in', 'FontSize', fntsz);\n");
fprintf(dat, "xlabel(['%c [' char(956) 'm]'], 'FontSize', fntsz); ylabel(['%c [' char(956) 'm]'], 'FontSize', fntsz);\n", zchr, xchr);
fprintf(dat, "global tle;\n");
fprintf(dat, "tle = title(' ', 'FontSize', fntsz+2);\n");
fprintf(dat, "pl_tih = text(%cend-(%cend-%cbeg)/30, %cbeg+(%cend-%cbeg)/30, ...\n", zchr, zchr, zchr, xchr, xchr, xchr);
fprintf(dat, "              sprintf('t = %%#4.2f T', pl_tp), 'Color', 'w', ...\n");
fprintf(dat, "              'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "              'VerticalAlignment', 'bottom', 'FontSize', fntsz);\n");
fprintf(dat, "hold on; box on; \n");
fprintf(dat, "\n");
fprintf(dat, "%% The initial plot\n");
fprintf(dat, "set(pl_img, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "for idx=1:nhd\n");
fprintf(dat, "   plot([hdbeg(idx), hdend(idx)], [hdpos(idx), hdpos(idx)], ...\n");
fprintf(dat, "        'w-', 'LineWidth', 1.0, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "end\n");
fprintf(dat, "for idx=1:nvd\n");
fprintf(dat, "   plot([vdpos(idx), vdpos(idx)], [vdbeg(idx), vdend(idx)], ...\n");
fprintf(dat, "        'w-', 'LineWidth', 1.0, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "end\n");
fprintf(dat, "ra = [0:2*pi/500:2*pi];\n");
fprintf(dat, "for idx=1:ncd\n");
fprintf(dat, "   plot(cos(ra)*cdrad(idx)+cdhp(idx), sin(ra)*cdrad(idx)+cdvp(idx), ...\n");
fprintf(dat, "        'w-', 'LineWidth', 1.0, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "end\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Crosshair \n");
fprintf(dat, "global chh chv ch%c ch%c;\n", xchr, zchr);
fprintf(dat, "ch%c = %cend; ch%c = %cend;\n", xchr, xchr, zchr, zchr);
fprintf(dat, "axes(axs);\n");
fprintf(dat, "chh = plot([%cbeg %cend], [ch%c ch%c], ...\n", zchr, zchr, xchr, xchr);
fprintf(dat, "           'LineStyle', '-.', 'LineWidth', 0.5, 'Color', [0.5 0.5 0.5], ...\n");
fprintf(dat, "           'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "chv = plot([ch%c ch%c], [%cbeg %cend], ...\n", zchr, zchr, xchr, xchr);
fprintf(dat, "           'LineStyle', '-.', 'LineWidth', 0.5, 'Color', [0.5 0.5 0.5], ...\n");
fprintf(dat, "           'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "\n");
fprintf(dat, "%% The colormap axis\n");
fprintf(dat, "global clb clb_img;\n");
fprintf(dat, "clb = axes;\n");
fprintf(dat, "rge = linspace(pl_min, pl_max, 64);\n");
fprintf(dat, "clb_img = imagesc([0 1], rge, [rge',rge'], [pl_min, pl_max]);\n");
fprintf(dat, "set(clb, 'Position', [cmpx axpy cmpw axh], 'Units', 'normalized', ...\n");
fprintf(dat, "         'YDir', 'normal', 'ButtonDownFcn', {@clbf_contour}, ...\n");
fprintf(dat, "         'XLim', [0, 1], 'YLim', [pl_min, pl_max], 'XTick', [], ...\n");
fprintf(dat, "         'YAxisLocation', 'right', 'TickDir', 'in', 'FontSize', fntsz);\n");
fprintf(dat, "set(clb_img, 'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "\n");
fprintf(dat, "%% The contour positioner\n");
fprintf(dat, "global cont_ah cont_minh cont_maxh cont_hdl cont_num; \n");
fprintf(dat, "lv = pl_min;\n");
fprintf(dat, "cont_ah = plot([0 1], [lv lv], ...\n");
fprintf(dat, "               'LineStyle', '-', 'LineWidth', 0.5, 'Color', [0.5 0.5 0.5], ...\n");
fprintf(dat, "               'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "cont_minh = plot([0 1], [pl_amin pl_amin], ...\n");
fprintf(dat, "                 'LineStyle', '-', 'LineWidth', 0.5, 'Color', [1.0 1.0 1.0], ...\n");
fprintf(dat, "                 'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "cont_maxh = plot([0 1], [pl_amax pl_amax], ...\n");
fprintf(dat, "                 'LineStyle', '-', 'LineWidth', 0.5, 'Color', [1.0 1.0 1.0], ...\n");
fprintf(dat, "                 'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "cont_num = 0;\n");
fprintf(dat, "\n");
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% The gui control objects \n");
fprintf(dat, "\n");
fprintf(dat, "uc_dummy = uicontrol('Style', 'text', 'String', 'M', 'Units', 'normalized', ...\n");
fprintf(dat, "                     'FontSize', fntsz, 'Visible', 'off');\n");
fprintf(dat, "dext = get(uc_dummy, 'Extent'); delete(uc_dummy);\n");
fprintf(dat, "uch0 = dext(4);\n");
fprintf(dat, "uch = (axh-37*busep)/11;\n");
fprintf(dat, "if uch0>uch uch0 = uch; end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotation: %c-coordinate \n", zchr);
fprintf(dat, "uch = uch0; ucpy = axpy; ucpx = bupx;\n");
fprintf(dat, "uc_zl = uicontrol('Style', 'text', 'String', '%c/um', 'Units', 'normalized', ...\n", zchr);
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "p = get(uc_zl, 'Extent'); ucw = p(3);\n");
fprintf(dat, "set(uc_zl, 'Position', [ucpx ucpy ucw uch]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Input field: %c-coordinate \n", zchr);
fprintf(dat, "global uc_zf;\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = buw-busep-ucw-busep-(buw-2*busep)/3;\n");
fprintf(dat, "uc_zf = uicontrol('Style', 'edit', 'Callback', {@cbf_zfld}, ...\n");
fprintf(dat, "                  'String', num2str(ch%c), 'Min', 1, 'Max', 1, ...\n", zchr);
fprintf(dat, "                  'TooltipString', 'View & set z-position', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_abg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: vertical cross section\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_vc = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_crosssec, 'v'}, ...\n");
fprintf(dat, "                  'String', '|', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Vertical field section', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotation: %c-coordinate \n", xchr);
fprintf(dat, "ucpy = ucpy+uch+busep; ucpx = bupx;\n");
fprintf(dat, "uc_xl = uicontrol('Style', 'text', 'String', '%c/um', 'Units', 'normalized', ...\n", xchr);
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "p = get(uc_xl, 'Extent'); ucw = p(3);\n");
fprintf(dat, "set(uc_xl, 'Position', [ucpx ucpy ucw uch]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Input field: %c-coordinate \n", xchr);
fprintf(dat, "global uc_xf;\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = buw-busep-ucw-busep-(buw-2*busep)/3;\n");
fprintf(dat, "uc_xf = uicontrol('Style', 'edit', 'Callback', {@cbf_xfld}, ...\n");
fprintf(dat, "                  'String', num2str(ch%c), 'Min', 1, 'Max', 1, ...\n", xchr);
fprintf(dat, "                  'TooltipString', 'View & set x-position', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_abg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: horizontal cross section\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_hc = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_crosssec, 'h'}, ...\n");
fprintf(dat, "                  'String', '---', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Horizontal field section', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: w-profile \n");
fprintf(dat, "ucw = (buw-2*busep)/3; ucpx = bupx; ucpy = ucpy+uch+10*busep;\n");
fprintf(dat, "uc_w  = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_wpl}, ...\n");
fprintf(dat, "                  'String', 'w', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Energy density of the optical field', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: refractive index profile\n");
fprintf(dat, "ucpx = bupx+2*ucw+2*busep;\n");
fprintf(dat, "uc_n  = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_npl}, ...\n");
fprintf(dat, "                  'String', 'n', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Refractive index profile', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Sx-profile \n");
fprintf(dat, "ucpx = bupx; ucpy = ucpy+uch+busep; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_sx = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_sxpl}, ...\n");
fprintf(dat, "                  'String', 'Sx', 'FontSize', fntsz, ...\n");
fprintf(dat, "            'TooltipString', 'Time-averaged Poynting vector, x-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Sy-profile \n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_sy = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_sypl}, ...\n");
fprintf(dat, "                  'String', 'Sy', 'FontSize', fntsz, ...\n");
fprintf(dat, "            'TooltipString', 'Time-averaged Poynting vector, y-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Sz-profile \n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_sz = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_szpl}, ...\n");
fprintf(dat, "                  'String', 'Sz', 'FontSize', fntsz, ...\n");
fprintf(dat, "            'TooltipString', 'Time-averaged Poynting vector, z-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: field-profile \n");
fprintf(dat, "ucpx = bupx; ucpy = ucpy+uch+busep; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_f  = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_fldpl}, ...\n");
fprintf(dat, "                  'String', 'F', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'The physical field, time snapshot', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: profile of the absolute field value\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_fm = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_fldmpl}, ...\n");
fprintf(dat, "                  'String', '|F|', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Absolute value of the complex field', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: profile of the absolute field square \n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_fq = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_fldqpl}, ...\n");
fprintf(dat, "                  'String', '|F|^2', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Absolute square of the complex field', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotations, slider time position\n");
fprintf(dat, "ucpy = ucpy+uch+busep; uch = uch0; ucw = buw/3; ucpx = bupx;\n");
fprintf(dat, "uc_tll = uicontrol('Style', 'text', 'String', 't=0 ...', ...\n");
fprintf(dat, "                   'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                   'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                   'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "ucpx = ucpx+ucw;\n");
fprintf(dat, "uc_tlc = uicontrol('Style', 'text', 'String', 'T/2', ...\n");
fprintf(dat, "                   'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                   'FontSize', fntsz, 'HorizontalAlignment', 'center', ...\n");
fprintf(dat, "                   'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "ucpx = ucpx+ucw;\n");
fprintf(dat, "uc_tlc = uicontrol('Style', 'text', 'String', '... T ', ...\n");
fprintf(dat, "                   'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                   'FontSize', fntsz, 'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "                   'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Slider: field snapshot time \n");
fprintf(dat, "global uc_tp;\n");
fprintf(dat, "ucpy = ucpy+uch+busep; uch = uch0*3/4; ucw = buw; ucpx = bupx;\n");
fprintf(dat, "uc_tp = uicontrol('Style', 'slider', 'Callback', {@cbf_timep}, ...\n");
fprintf(dat, "                  'Max', 1.0, 'Min', 0.0, 'Value', 0.0, ...\n");
fprintf(dat, "      'TooltipString', 'Physical field, snapshot time within one period T', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: play reverse animation\n");
fprintf(dat, "ucpy = ucpy+uch+busep; uch = uch0; ucpx = bupx; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_anbw = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_play, 'b'}, ...\n");
fprintf(dat, "                    'String', '<<', 'FontSize', fntsz, ...\n");
fprintf(dat, "                    'TooltipString', 'Backward wave propagation', ...\n");
fprintf(dat, "                    'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                    'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Input field: number of time periods shown in the animation\n");
fprintf(dat, "global uc_annp; ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_annp = uicontrol('Style', 'edit', 'Callback', {@cbf_annp}, ...\n");
fprintf(dat, "                    'String', num2str(an_nump), 'Min', 1, 'Max', 1, ...\n");
fprintf(dat, "                    'TooltipString', 'Runtime/T of the animation', ...\n");
fprintf(dat, "                    'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                    'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                    'BackgroundColor', col_abg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: play field animation\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_anfw = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_play, 'f'}, ...\n");
fprintf(dat, "                    'String', '>>', 'FontSize', fntsz, ...\n");
fprintf(dat, "                    'TooltipString', 'Forward wave propagation', ...\n");
fprintf(dat, "                    'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                    'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Hx-component\n");
fprintf(dat, "ucpy = ucpy+uch+10*busep; uch = uch0; ucpx = bupx; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_hx = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_hxpl}, ...\n");
fprintf(dat, "                  'String', 'Hx', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Magnetic field, x-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Hy-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_hy = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_hypl}, ...\n");
fprintf(dat, "                  'String', 'Hy', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Magnetic field, y-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Hz-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_hz = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_hzpl}, ...\n");
fprintf(dat, "                  'String', 'Hz', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Magnetic field, z-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Ex-component\n");
fprintf(dat, "ucpx = bupx; ucpy = ucpy+uch+busep; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_ex = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_expl}, ...\n");
fprintf(dat, "                  'String', 'Ex', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Electric field, x-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Ey-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_ey = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_eypl}, ...\n");
fprintf(dat, "                  'String', 'Ey', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Electric field, y-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Ey-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_ez = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_ezpl}, ...\n");
fprintf(dat, "                  'String', 'Ez', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Electric field, z-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotation: polarization and wavelength\n");
fprintf(dat, "uch = uch0; ucpy = axpy+axh-uch; ucpx = bupx;\n");
fprintf(dat, "uc_an = uicontrol('Style', 'text', 'Units', 'normalized', ...\n");
fprintf(dat, "                  'String', [pol ', ' num2str(lambda) ' um'], ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "p = get(uc_an, 'Extent'); ucw = p(3);\n");
fprintf(dat, "set(uc_an, 'Position', [ucpx ucpy ucw uch]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% Callback functions for the gui objects \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: input z-coordinate \n");
fprintf(dat, "function cbf_zfld(obj, eventdata)\n");
fprintf(dat, "global uc_zf ch%c chv %cbeg %cend;\n", zchr, zchr, zchr);
fprintf(dat, "ch%c = str2double(get(uc_zf, 'String'));\n", zchr);
fprintf(dat, "if isnan(ch%c) ch%c = 0.5*(%cbeg+%cend); end\n", zchr, zchr, zchr, zchr);
fprintf(dat, "if ch%c<%cbeg ch%c = %cbeg; end\n", zchr, zchr, zchr, zchr);
fprintf(dat, "if ch%c>%cend ch%c = %cend; end\n", zchr, zchr, zchr, zchr);
fprintf(dat, "set(chv, 'XData', [ch%c ch%c]); \n", zchr, zchr);
fprintf(dat, "set(uc_zf, 'String', num2str(ch%c)); \n", zchr);
fprintf(dat, "\n");
fprintf(dat, "%% Callback: input x-coordinate \n");
fprintf(dat, "function cbf_xfld(obj, eventdata)\n");
fprintf(dat, "global uc_xf ch%c chh %cbeg %cend;\n", xchr, xchr, xchr);
fprintf(dat, "ch%c = str2double(get(uc_xf, 'String'));\n", xchr);
fprintf(dat, "if isnan(ch%c) ch%c = 0.5*(%cbeg+%cend); end\n", xchr, xchr, xchr, xchr);
fprintf(dat, "if ch%c<%cbeg ch%c = %cbeg; end\n", xchr, xchr, xchr, xchr);
fprintf(dat, "if ch%c>%cend ch%c = %cend; end\n", xchr, xchr, xchr, xchr);
fprintf(dat, "set(chh, 'YData', [ch%c ch%c]); \n", xchr, xchr);
fprintf(dat, "set(uc_xf, 'String', num2str(ch%c)); \n", xchr);
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Sx \n");
fprintf(dat, "function cbf_sxpl(obj, eventdata)\n");
fprintf(dat, "global Sx max_S pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Sx; pl_fldim = Sx*0; pl_cp = 'S_x';\n");
fprintf(dat, "pl_min = -max_S; pl_max = max_S;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Sy \n");
fprintf(dat, "function cbf_sypl(obj, eventdata)\n");
fprintf(dat, "global Sy max_S pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Sy; pl_fldim = Sy*0; pl_cp = 'S_y';\n");
fprintf(dat, "pl_min = -max_S; pl_max = max_S;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Sz \n");
fprintf(dat, "function cbf_szpl(obj, eventdata)\n");
fprintf(dat, "global Sz max_S pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Sz; pl_fldim = Sz*0; pl_cp = 'S_z';\n");
fprintf(dat, "pl_min = -max_S; pl_max = max_S;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select w \n");
fprintf(dat, "function cbf_wpl(obj, eventdata)\n");
fprintf(dat, "global W max_W pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = W; pl_fldim = W*0; pl_cp = 'w'; pl_min = 0.0; pl_max = max_W;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select physical field \n");
fprintf(dat, "function cbf_fldpl(obj, eventdata)\n");
fprintf(dat, "global pl_afo;\n");
fprintf(dat, "pl_afo = 'f'; \n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select absolute field value\n");
fprintf(dat, "function cbf_fldmpl(obj, eventdata)\n");
fprintf(dat, "global pl_afo;\n");
fprintf(dat, "pl_afo = 'm'; \n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select absolute field square \n");
fprintf(dat, "function cbf_fldqpl(obj, eventdata)\n");
fprintf(dat, "global pl_afo;\n");
fprintf(dat, "pl_afo = 'q'; \n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select adjust snapshot time \n");
fprintf(dat, "function cbf_timep(obj, eventdata)\n");
fprintf(dat, "global uc_tp pl_tp;\n");
fprintf(dat, "pl_tp = get(uc_tp, 'Value');\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select animation runtime \n");
fprintf(dat, "function cbf_annp(obj, eventdata)\n");
fprintf(dat, "global uc_annp an_nump;\n");
fprintf(dat, "np = str2double(get(uc_annp, 'String'));\n");
fprintf(dat, "if isnan(np) np = 10; end\n");
fprintf(dat, "if np<1  np = 1;  end\n");
fprintf(dat, "if np>99 np = 99; end\n");
fprintf(dat, "an_nump = round(np);\n");
fprintf(dat, "set(uc_annp, 'String', num2str(an_nump)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Hx \n");
fprintf(dat, "function cbf_hxpl(obj, eventdata)\n");
fprintf(dat, "global Hx_re Hx_im max_H pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Hx_re; pl_fldim = Hx_im; pl_cp = 'H_x';\n");
fprintf(dat, "pl_min = -max_H; pl_max = max_H;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Hy \n");
fprintf(dat, "function cbf_hypl(obj, eventdata)\n");
fprintf(dat, "global Hy_re Hy_im max_H pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Hy_re; pl_fldim = Hy_im; pl_cp = 'H_y';\n");
fprintf(dat, "pl_min = -max_H; pl_max = max_H;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Hz \n");
fprintf(dat, "function cbf_hzpl(obj, eventdata)\n");
fprintf(dat, "global Hz_re Hz_im max_H pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Hz_re; pl_fldim = Hz_im; pl_cp = 'H_z';\n");
fprintf(dat, "pl_min = -max_H; pl_max = max_H;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Ex \n");
fprintf(dat, "function cbf_expl(obj, eventdata)\n");
fprintf(dat, "global Ex_re Ex_im max_E pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Ex_re; pl_fldim = Ex_im; pl_cp = 'E_x';\n");
fprintf(dat, "pl_min = -max_E; pl_max = max_E;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Ey \n");
fprintf(dat, "function cbf_eypl(obj, eventdata)\n");
fprintf(dat, "global Ey_re Ey_im max_E pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Ey_re; pl_fldim = Ey_im; pl_cp = 'E_y';\n");
fprintf(dat, "pl_min = -max_E; pl_max = max_E;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Ez \n");
fprintf(dat, "function cbf_ezpl(obj, eventdata)\n");
fprintf(dat, "global Ez_re Ez_im max_E pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Ez_re; pl_fldim = Ez_im; pl_cp = 'E_z';\n");
fprintf(dat, "pl_min = -max_E; pl_max = max_E;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select n\n");
fprintf(dat, "function cbf_npl(obj, eventdata)\n");
fprintf(dat, "global n min_n max_n pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = n; pl_fldim = n*0.0; pl_cp = 'n'; pl_min = min_n; pl_max = max_n;\n");
fprintf(dat, "if pl_min>=pl_max pl_min = pl_min-0.5; pl_max = pl_min+1.0; end;\n");
fprintf(dat, "d = (pl_max-pl_min)/3; pl_min = pl_min-d; pl_max = pl_max+d;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% The actual plotting routines, partly callbacks for gui objects\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: crosshair positioning \n");
fprintf(dat, "function clbf_crossh(obj, eventdata)\n");
fprintf(dat, "global fgr axs chh chv ch%c ch%c;\n", xchr, zchr);
fprintf(dat, "mpos = get(0, 'PointerLocation'); \n");
fprintf(dat, "mx=mpos(1); my=mpos(2);\n");
fprintf(dat, "fpos = get(fgr, 'Position'); fx=fpos(1); fy=fpos(2); fw=fpos(3); fh=fpos(4); \n");
fprintf(dat, "rx = (mx-fx)/fw; ry = (my-fy)/fh;\n");
fprintf(dat, "apos = get(axs, 'Position'); ax=apos(1); ay=apos(2); aw=apos(3); ah=apos(4); \n");
fprintf(dat, "extx = get(axs, 'Xlim'); x0=extx(1); x1=extx(2);\n");
fprintf(dat, "exty = get(axs, 'Ylim'); y0=exty(1); y1=exty(2);\n");
fprintf(dat, "px = x0+(rx-ax)/aw*(x1-x0); py = y0+(ry-ay)/ah*(y1-y0);\n");
fprintf(dat, "ch%c = py; ch%c = px;\n", xchr, zchr);
fprintf(dat, "global %cbeg %cend %cbeg %cend;\n", xchr, xchr, zchr, zchr);
fprintf(dat, "if ch%c<%cbeg ch%c = %cbeg; end\n", xchr, xchr, xchr, xchr);
fprintf(dat, "if ch%c>%cend ch%c = %cend; end\n", xchr, xchr, xchr, xchr);
fprintf(dat, "if ch%c<%cbeg ch%c = %cbeg; end\n", zchr, zchr, zchr, zchr);
fprintf(dat, "if ch%c>%cend ch%c = %cend; end\n", zchr, zchr, zchr, zchr);
fprintf(dat, "set(chh, 'YData', [ch%c ch%c]); set(chv, 'XData', [ch%c ch%c]); \n", xchr, xchr, zchr, zchr);
fprintf(dat, "global uc_xf uc_zf;\n");
fprintf(dat, "set(uc_xf, 'String', num2str(ch%c)); set(uc_zf, 'String', num2str(ch%c)); \n", xchr, zchr);
fprintf(dat, "\n");
fprintf(dat, "%% Callback: place contour line\n");
fprintf(dat, "function clbf_contour(obj, eventdata)\n");
fprintf(dat, "global fgr axs clb cont_ah cont_num cont_hdl;\n");
fprintf(dat, "for k=1:cont_num delete(cont_hdl(k)); end\n");
fprintf(dat, "mpos = get(0, 'PointerLocation'); my=mpos(2);\n");
fprintf(dat, "fpos = get(fgr, 'Position'); fy=fpos(2); fh=fpos(4); \n");
fprintf(dat, "ry = (my-fy)/fh;\n");
fprintf(dat, "apos = get(clb, 'Position'); ay=apos(2); ah=apos(4); \n");
fprintf(dat, "exty = get(clb, 'Ylim'); y0=exty(1); y1=exty(2);\n");
fprintf(dat, "lv = y0+(ry-ay)/ah*(y1-y0); \n");
fprintf(dat, "global pl_amin pl_amax;\n");
fprintf(dat, "if lv<pl_amin lv = pl_amin; end\n");
fprintf(dat, "if lv>pl_amax lv = pl_amax; end\n");
fprintf(dat, "set(cont_ah, 'YData', [lv lv]); \n");
fprintf(dat, "axes(axs); hold on;\n");
fprintf(dat, "global %c %c pl_fld;\n", zchr, xchr);
fprintf(dat, "cm = contourc(%c, %c, pl_fld, [lv, lv]); cl = size(cm, 2);\n", zchr, xchr);
fprintf(dat, "k = 1; cont_hdl = ones(1, cl); cont_num = 0;\n");
fprintf(dat, "while(k < cl)\n");
fprintf(dat, "      np = cm(2, k); pz = cm(1, k+1:k+np); px = cm(2, k+1:k+np);\n");
fprintf(dat, "      cont_num = cont_num+1;\n");
fprintf(dat, "      cont_hdl(cont_num) = line(pz, px, ...\n");
fprintf(dat, "                          'Color', 'b', 'LineStyle', '-', 'LineWidth', 0.5, ...\n");
fprintf(dat, "                          'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "      k = k+np+1;\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: update the field plot \n");
fprintf(dat, "function updateplot\n");
fprintf(dat, "global fgr axs tle clb clb_img;\n");
fprintf(dat, "global pl_fld pl_fldre pl_fldim pl_tp pl_min pl_max pl_amin pl_amax;\n");
fprintf(dat, "global pl_cp pl_img pl_afo pl_tih;\n");
fprintf(dat, "global %c %c;\n", xchr, zchr);
fprintf(dat, "global cont_ah cont_minh cont_maxh cont_num cont_hdl;\n");
fprintf(dat, "global mag_blue lin_blue;\n");
fprintf(dat, "for k=1:cont_num delete(cont_hdl(k)); end\n");
fprintf(dat, "cont_num = 0;\n");
fprintf(dat, "switch pl_cp\n");
fprintf(dat, "case {'E_x', 'E_y', 'E_z', 'H_x', 'H_y', 'H_z'} \n");
fprintf(dat, "   if pl_afo == 'm' \n");
fprintf(dat, "      title = ['| ' pl_cp '(%c, %c) |']; \n", xchr, zchr);
fprintf(dat, "      pl_fld = sqrt(pl_fldre.*pl_fldre+pl_fldim.*pl_fldim);\n");
fprintf(dat, "      clm = lin_blue; minlv = 0.0; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   elseif pl_afo == 'q'\n");
fprintf(dat, "      title = ['| ' pl_cp '(%c, %c) |^2']; \n", xchr, zchr);
fprintf(dat, "      pl_fld = pl_fldre.*pl_fldre+pl_fldim.*pl_fldim;\n");
fprintf(dat, "      clm = lin_blue; minlv = 0.0; maxlv = pl_max*pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   else\n");
fprintf(dat, "      title = ['Re ' pl_cp '(%c, %c) e^{i2' char(960) 't/T}'];\n", xchr, zchr);
fprintf(dat, "      pl_fld = pl_fldre*cos(2*pi*pl_tp)-pl_fldim*sin(2*pi*pl_tp);\n");
fprintf(dat, "      clm = mag_blue; minlv = -pl_max; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'on', 'String', sprintf('t = %%#4.2f T', pl_tp));\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'S_x', 'S_y', 'S_z'}\n");
fprintf(dat, "   if pl_afo == 'm' || pl_afo == 'q'\n");
fprintf(dat, "      title = ['| ' pl_cp '(%c, %c) |']; \n", xchr, zchr);
fprintf(dat, "      pl_fld = sqrt(pl_fldre.*pl_fldre);\n");
fprintf(dat, "      clm = lin_blue; minlv = 0.0; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   else\n");
fprintf(dat, "      title = [pl_cp '(%c, %c)'];\n", xchr, zchr);
fprintf(dat, "      pl_fld = pl_fldre;\n");
fprintf(dat, "      clm = mag_blue; minlv = -pl_max; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'w'}\n");
fprintf(dat, "   title = 'w(%c, %c)';\n", xchr, zchr);
fprintf(dat, "   pl_fld = pl_fldre;\n");
fprintf(dat, "   clm = lin_blue; minlv = 0.0; maxlv = pl_max;\n");
fprintf(dat, "   set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "case {'n'}\n");
fprintf(dat, "   title = 'n(%c, %c)';\n", xchr, zchr);
fprintf(dat, "   pl_fld = pl_fldre;\n");
fprintf(dat, "   clm = lin_blue; minlv = pl_min; maxlv = pl_max;\n");
fprintf(dat, "   set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "end\n");
fprintf(dat, "if minlv == maxlv\n");
fprintf(dat, "   minlv = minlv-1.0\n");
fprintf(dat, "   maxlv = minlv+2.0\n");
fprintf(dat, "end\n");
fprintf(dat, "pl_amin = min(min(pl_fld)); pl_amax = max(max(pl_fld));\n");
fprintf(dat, "set(fgr, 'Colormap', clm);\n");
fprintf(dat, "set(clb_img, 'CData', linspace(minlv, maxlv, 64)', ...\n");
fprintf(dat, "             'YData', linspace(minlv, maxlv, 64)');\n");
fprintf(dat, "set(clb, 'YLim', [minlv, maxlv], 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "set(cont_ah, 'YData', [minlv minlv]); \n");
fprintf(dat, "set(cont_minh, 'YData', [pl_amin pl_amin]); \n");
fprintf(dat, "set(cont_maxh, 'YData', [pl_amax pl_amax]); \n");
fprintf(dat, "set(tle, 'String', title);\n");
fprintf(dat, "set(axs, 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "set(pl_img, 'CData', pl_fld);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: field animation \n");
fprintf(dat, "function cbf_play(obj, eventdata, dr) \n");
fprintf(dat, "global fgr axs tle clb clb_img uc_tp mag_blue lin_blue;\n");
fprintf(dat, "global pl_fld pl_fldre pl_fldim pl_tp pl_min pl_max %c %c;\n", xchr, zchr);
fprintf(dat, "global pl_cp pl_img pl_afo pl_tih an_nump;\n");
fprintf(dat, "global cont_ah cont_num cont_hdl;\n");
fprintf(dat, "for k=1:cont_num delete(cont_hdl(k)); end\n");
fprintf(dat, "cont_num = 0;\n");
fprintf(dat, "switch pl_cp\n");
fprintf(dat, "case {'E_x', 'E_y', 'E_z', 'H_x', 'H_y', 'H_z'} \n");
fprintf(dat, "   if pl_afo == 'm' || pl_afo =='q' return; end\n");
fprintf(dat, "case {'S_x', 'S_y', 'S_z', 'w', 'n'}\n");
fprintf(dat, "   return;\n");
fprintf(dat, "end\n");
fprintf(dat, "clm = mag_blue;\n");
fprintf(dat, "minlv = -pl_max; maxlv = pl_max;\n");
fprintf(dat, "set(fgr, 'Colormap', clm);\n");
fprintf(dat, "set(clb_img, 'CData', linspace(minlv, maxlv, 64)', ...\n");
fprintf(dat, "             'YData', linspace(minlv, maxlv, 64)');\n");
fprintf(dat, "set(clb, 'YLim', [minlv, maxlv], 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "set(cont_ah, 'YData', [minlv minlv]); \n");
fprintf(dat, "set(axs, 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "title = ['Re ' pl_cp '(%c, %c) e^{i2' char(960) 't/T}'];\n", xchr, zchr);
fprintf(dat, "set(tle, 'String', title);\n");
fprintf(dat, "set(pl_tih, 'Visible', 'on', 'String', sprintf('t = %%#4.2f T', pl_tp));\n");
fprintf(dat, "axes(axs);\n");
fprintf(dat, "for li=1:an_nump\n");
fprintf(dat, "   for ti=0:29\n");
fprintf(dat, "      tic;\n");
fprintf(dat, "      itim = toc;\n");
fprintf(dat, "      if dr == 'b'\n");
fprintf(dat, "         t = pl_tp-ti/30; if t<0 t = 1+t; end\n");
fprintf(dat, "      else \n");
fprintf(dat, "         t = pl_tp+ti/30; if t>1 t = t-1; end\n");
fprintf(dat, "      end\n");
fprintf(dat, "      pl_fld = pl_fldre*cos(2*pi*t)-pl_fldim*sin(2*pi*t);\n");
fprintf(dat, "      set(pl_img, 'CData', pl_fld);\n");
fprintf(dat, "      set(pl_tih, 'String', sprintf('t = %%#4.2f T', t));\n");
fprintf(dat, "      set(uc_tp, 'Value', t);\n");
fprintf(dat, "      drawnow;\n");
fprintf(dat, "      etim = toc;\n");
fprintf(dat, "      while etim-itim < 0.0666\n");
fprintf(dat, "         etim = toc;\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "set(uc_tp, 'Value', pl_tp);\n");
fprintf(dat, "title = ['Re ' pl_cp '(%c, %c) e^{i2' char(960) 't/T}'];\n", xchr, zchr);
fprintf(dat, "set(tle, 'String', title);\n");
fprintf(dat, "pl_fld = pl_fldre*cos(2*pi*(pl_tp+t))-pl_fldim*sin(2*pi*(pl_tp+t));\n");
fprintf(dat, "set(pl_img, 'CData', pl_fld);\n");
fprintf(dat, "set(pl_tih, 'String', sprintf('t = %%#4.2f T', pl_tp));\n");
fprintf(dat, "clear msp;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: cross section figure\n");
fprintf(dat, "function cbf_crosssec(obj, eventdata, dr)\n");
fprintf(dat, "global csf_xaw csf_zaw;\n");
fprintf(dat, "global col_fbg col_txt col_abg col_atx fntsz fgr;\n");
fprintf(dat, "global %c %c %cbeg %cend %cbeg %cend num%c num%c;\n", xchr, zchr, xchr, xchr, zchr, zchr, xchr, zchr);
fprintf(dat, "global ch%c ch%c pl_fld pl_cp pl_afo pl_tp pl_min pl_max; \n", xchr, zchr);
fprintf(dat, "csf = figure('NumberTitle', 'off', 'Name', '%s, field section', ...\n", name);
fprintf(dat, "             'Color', col_fbg, 'Units', 'pixels', 'Menubar', 'none');\n");
fprintf(dat, "if dr == 'h' axw = csf_zaw; else axw = csf_xaw; end\n");
fprintf(dat, "axh = 200; fgw = axw+60+20; fgh = axh+40+20;\n");
fprintf(dat, "axw = axw/fgw; axh = axh/fgh; axpx = 60/fgw; axpy = 40/fgh;\n");
fprintf(dat, "scrsz = get(0, 'ScreenSize');\n");
fprintf(dat, "if fgw > scrsz(3)*0.6 \n");
fprintf(dat, "   s = scrsz(3)*0.6/fgw; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "if fgh > scrsz(4)*0.6 \n");
fprintf(dat, "   s = scrsz(4)*0.6/fgh; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "fgpx = round(scrsz(1)+scrsz(3)*0.1); fgpy = round(scrsz(2)+scrsz(4)*0.90)-fgh;\n");
fprintf(dat, "p = get(fgr, 'Position');\n");
fprintf(dat, "if p(1)+p(3)+fgw+10 <= scrsz(1)+scrsz(3) \n");
fprintf(dat, "   fgpx = p(1)+p(3)+10; fgpy = p(2);\n");
fprintf(dat, "end\n");
fprintf(dat, "set(csf, 'Position', [fgpx fgpy fgw fgh]);\n");
fprintf(dat, "csf_axs = axes;\n");
fprintf(dat, "axes(csf_axs); hold on; box on;\n");
fprintf(dat, "set(csf_axs, 'Position', [axpx axpy axw axh], 'Units', 'normalized', ...\n");
fprintf(dat, "             'YDir', 'normal', 'TickDir', 'out', 'FontSize', fntsz);\n");
fprintf(dat, "if dr == 'h'\n");
fprintf(dat, "   xlabel(['%c [' char(956) 'm]'], 'FontSize', fntsz);\n", zchr);
fprintf(dat, "   cbeg = %cbeg; cend = %cend; cod = %c;\n", zchr, zchr, zchr);
// fprintf(dat, "   poss = sprintf('%c = %%.5g \\\\mum', ch%c);\n", xchr, xchr);
fprintf(dat, "   poss = sprintf(['%c = %%.5g ' char(956) 'm'], ch%c);\n", xchr, xchr);
fprintf(dat, "else \n");
fprintf(dat, "   xlabel(['%c [' char(956) 'm]'], 'FontSize', fntsz);\n", xchr);
fprintf(dat, "   cbeg = %cbeg; cend = %cend; cod = %c;\n", xchr, xchr, xchr);
// fprintf(dat, "   poss = sprintf('%c = %%.5g \\\\mum', ch%c);\n", zchr, zchr);
fprintf(dat, "   poss = sprintf(['%c = %%.5g ' char(956) 'm'], ch%c);\n", zchr, zchr);
fprintf(dat, "end\n");
fprintf(dat, "set(csf_axs, 'XLim', [cbeg, cend]);\n");
fprintf(dat, "coc = linspace(cbeg, cend, 1000);\n");
fprintf(dat, "if dr == 'h'\n");
fprintf(dat, "   fldicoc = interp2(%c, %c, pl_fld, coc, linspace(ch%c, ch%c, 1000), 'linear');\n", zchr, xchr, xchr, xchr);
fprintf(dat, "   fldicod = interp2(%c, %c, pl_fld, cod, linspace(ch%c, ch%c, length(%c)));\n", zchr, xchr, xchr, xchr, zchr);
fprintf(dat, "else\n");
fprintf(dat, "   fldicoc = interp2(%c, %c, pl_fld, linspace(ch%c, ch%c, 1000), coc, 'linear');\n", zchr, xchr, zchr, zchr);
fprintf(dat, "   fldicod = interp2(%c, %c, pl_fld, linspace(ch%c, ch%c, length(%c)), cod);\n", zchr, xchr, zchr, zchr, xchr);
fprintf(dat, "end\n");
fprintf(dat, "switch pl_cp\n");
fprintf(dat, "case {'E_x', 'E_y', 'E_z', 'H_x', 'H_y', 'H_z'} \n");
fprintf(dat, "   if pl_afo == 'm' \n");
fprintf(dat, "      ylabel(['| ' pl_cp ' |'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = 0.0; fend = pl_max;\n");
fprintf(dat, "   elseif pl_afo == 'q'\n");
fprintf(dat, "      ylabel(['| ' pl_cp ' |^2'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = 0.0; fend = pl_max*pl_max;\n");
fprintf(dat, "   else\n");
fprintf(dat, "      ylabel(['Re ' pl_cp ' e^{i2' char(960) 't/T}'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = -pl_max; fend = pl_max;\n");
fprintf(dat, "      poss = [poss, sprintf(',  t = %%#4.2f T', pl_tp)];\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'S_x', 'S_y', 'S_z'}\n");
fprintf(dat, "   if pl_afo == 'm' || pl_afo == 'q'\n");
fprintf(dat, "      ylabel(['| ' pl_cp ' |'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = 0.0; fend = pl_max;\n");
fprintf(dat, "   else\n");
fprintf(dat, "      ylabel(pl_cp, 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = -pl_max; fend = pl_max;\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'w'}\n");
fprintf(dat, "   ylabel('w', 'FontSize', fntsz);\n");
fprintf(dat, "   fbeg = 0.0; fend = pl_max;\n");
fprintf(dat, "case {'n'}\n");
fprintf(dat, "   ylabel('n', 'FontSize', fntsz);\n");
fprintf(dat, "   fbeg = pl_min; fend = pl_max;\n");
fprintf(dat, "end\n");
fprintf(dat, "global n nvd vdbeg vdend vdpos nhd hdbeg hdend hdpos min_n max_n;\n");
fprintf(dat, "if dr == 'h'\n");
fprintf(dat, "   cb = ones(nvd); ncb = 0;\n");
fprintf(dat, "   for d=1:nvd\n");
fprintf(dat, "      if vdbeg(d) <= ch%c && ch%c < vdend(d)\n", xchr, xchr);
fprintf(dat, "         ncb = ncb+1; cb(ncb) = vdpos(d);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for d=1:ncb+1\n");
fprintf(dat, "      if d==1     b0=cbeg; else b0=cb(d-1); end\n");
fprintf(dat, "      if d==ncb+1 b1=cend; else b1=cb(d);   end\n");
fprintf(dat, "      rip = 0.5*(b1+b0);\n");
fprintf(dat, "      if %cbeg <= rip && rip <= %cend\n", zchr, zchr);
fprintf(dat, "         ri = interp2(%c, %c, n, 0.5*(b1+b0), ch%c);\n", zchr, xchr, xchr);
fprintf(dat, "         pc = 1.0-(ri-min_n)/(max_n-min_n)*0.2;\n");
fprintf(dat, "         patch([b0 b0 b1 b1], [fend fbeg fbeg fend], [pc pc pc], ...\n");
fprintf(dat, "              'LineStyle', 'none');\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "else\n");
fprintf(dat, "   cb = ones(nhd); ncb = 0;\n");
fprintf(dat, "   for d=1:nhd\n");
fprintf(dat, "      if hdbeg(d) <= ch%c && ch%c < hdend(d)\n", zchr, zchr);
fprintf(dat, "         ncb = ncb+1; cb(ncb) = hdpos(d);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for d=1:ncb+1\n");
fprintf(dat, "      if d==1     b0=cbeg; else b0=cb(d-1); end\n");
fprintf(dat, "      if d==ncb+1 b1=cend; else b1=cb(d);   end\n");
fprintf(dat, "      rip = 0.5*(b1+b0);\n");
fprintf(dat, "      if %cbeg <= rip && rip <= %cend\n", xchr, xchr);
fprintf(dat, "         ri = interp2(%c, %c, n, ch%c, 0.5*(b1+b0));\n", zchr, xchr, zchr);
fprintf(dat, "         pc = 1.0-(ri-min_n)/(max_n-min_n)*0.2;\n");
fprintf(dat, "         patch([b0 b0 b1 b1], [fend fbeg fbeg fend], [pc pc pc], ...\n");
fprintf(dat, "              'LineStyle', 'none');\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for d=1:ncb\n");
fprintf(dat, "   plot([cb(d), cb(d)], [fbeg fend], 'k-', 'LineWidth', 1);\n");
fprintf(dat, "end\n");
fprintf(dat, "text(cend-(cend-cbeg)/30, fend-(fend-fbeg)/30, ...\n");
fprintf(dat, "     poss, 'Color', 'k', 'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "     'VerticalAlignment', 'top', 'FontSize', fntsz);\n");
fprintf(dat, "set(csf_axs, 'YLim', [fbeg, fend]);\n");
fprintf(dat, "plot(coc, fldicoc, 'b-', 'LineWidth', 2); \n");
fprintf(dat, "ms = 2; if length(cod) <= 60 ms = 4; end\n");
fprintf(dat, "plot(cod, fldicod, 'ko', 'MarkerSize', ms, 'LineWidth', 2); \n");
fprintf(dat, "plot([cbeg cbeg cend cend cbeg], [fend fbeg fbeg fend fend], 'k-');\n");
fprintf(dat, "\n");
}

/* 3D waveguide cross section: geometry */
void mlout_gengeoxy(FILE *dat, SegWgCrs s, Rect disp)
{
	fprintf(dat, "%% Geometry\n");
	fprintf(dat, "\n");

	Dvector xpos;
	Dvector beg;
	Dvector end;
	Dvector pos;
	int x1, x2, p;
	double x1p, x2p, y0p, y1p, r; 
	int nd;

	fprintf(dat, "%% Vertical discontinuities\n");
	for(int j=0; j<=s.ny; ++j)
	{
		x1 = 0; x2 = 0; xpos.strip();
		while(x1 <= s(j).nx || x2 <= s(j+1).nx)
		{
			if(x1<=s(j).nx) x1p = s(j).hx(x1); 
			else x1p=1.0e+300;
			if(x2<=s(j+1).nx) x2p = s(j+1).hx(x2); 
			else x2p=1.0e+300;
			if(x1p < x2p)
			{
				xpos.append(x1p);
				++x1;
			}
			else
			{
				xpos.append(x2p);
				++x2;
			}
		}
		p = xpos.nel;
		for(int i=-1; i<=p-1; ++i)
		{
			if(i==-1) x1p = disp.x0; 
			else      x1p = xpos(i);
			if(i==p-1) x2p = disp.x1;
			else       x2p = xpos(i+1);
			if(fabs(x1p-x2p) >= 1.0e-10)
			{
				r = 0.5*(x1p+x2p);
				if(fabs(s(j).eps(r)-s(j+1).eps(r))
				   >1.0e-10)
				{
					beg.append(x1p);
					end.append(x2p);
					pos.append(s.hy(j));
				} 
			}
		}
	}
	nd = pos.nel;
	fprintf(dat, "nvd = %d;\n", nd);
	fprintf(dat, "vdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "vdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	pos.strip();
	beg.strip();
	end.strip();

	fprintf(dat, "%% Horizontal discontinuities\n");
	for(int j=0; j<=s.ny+1; ++j)
	{
		if(j == 0) y0p = disp.y0; 
		else       y0p = s.hy(j-1);
		if(j == s.ny+1) y1p = disp.y1; 
		else       	y1p = s.hy(j);
		for(int i=0; i<=s(j).nx; ++i)
		{
			if(fabs(s(j).eps(i)-s(j).eps(i+1))>1.0e-10)
			{
				pos.append(s(j).hx(i));
				beg.append(y0p);
				end.append(y1p);
			}
		}
	}
	nd = pos.nel;
	fprintf(dat, "nhd = %d;\n", nd);
	fprintf(dat, "hdbeg = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", beg(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdend = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", end(i));
	fprintf(dat, "];\n");
	fprintf(dat, "hdpos = [ ");
	for(int i=0; i<=nd-1; ++i) fprintf(dat, "%g ", pos(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	fprintf(dat, "%% Circular discontinuities\n");
	fprintf(dat, "ncd = %d;\n", 0);
	fprintf(dat, "cdrad = [ ");
	fprintf(dat, "];\n");
	fprintf(dat, "cdvp  = [ ");
	fprintf(dat, "];\n");
	fprintf(dat, "cdhp  = [ ");
	fprintf(dat, "];\n");
	fprintf(dat, "\n");

	return;
}

/* 3D waveguide cross section: refractive index profile patches, data */
void mlout_gengeorefindpatches(FILE *dat, SegWgCrs s) 
{
        fprintf(dat, "%% Refractive index patches\n");
	WgCrs c = s.wgcrs();
	double faraway = (c.hx(c.nx)-c.hx(0)+c.hy(c.ny)-c.hy(0))*10.0;
	int nump = (c.nx+2)*(c.ny+2);
	Dvector vx0(nump), vx1(nump), vy0(nump), vy1(nump), vn(nump);
	nump = 0;
	double x0, x1, y0, y1;
	for(int xi=0; xi<=c.nx+1; ++xi)
	{
		if(xi == 0) x0 = -faraway;
		else        x0 = c.hx(xi-1);
		if(xi == c.nx+1) x1 = faraway;
		else             x1 = c.hx(xi);
		for(int yi=0; yi<=c.ny+1; ++yi)
		{
			if(yi == 0) y0 = -faraway;
			else        y0 = c.hy(yi-1);
			if(yi == c.ny+1) y1 = faraway;
			else             y1 = c.hy(yi);
			vx0(nump) = x0;
			vx1(nump) = x1;
			vy0(nump) = y0;
			vy1(nump) = y1;
			vn(nump) = c.n(xi, yi); 
			++nump;
		}
	}
	fprintf(dat, "numnp = %d;\n", nump);
	fprintf(dat, "npx0 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vx0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npx1 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vx1(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npy0 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vy0(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npy1 = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vy1(i));
	fprintf(dat, "];\n");
	fprintf(dat, "npn = [ ");
	for(int i=0; i<=nump-1; ++i) fprintf(dat, "%g ", vn(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
}

/* output window, x-y-plane */
void mlout_windowxy(FILE *dat, double xbeg, double xend,
                               double ybeg, double yend)
{
	fprintf(dat, "%% Output window\n");
	fprintf(dat, "xbeg = %g;\n", xbeg);
	fprintf(dat, "xend = %g;\n", xend);
	fprintf(dat, "ybeg = %g;\n", ybeg);
	fprintf(dat, "yend = %g;\n", yend);
	fprintf(dat, "\n");
	return;
}

/* 2D waveguide cross ssection, plot commands */
void mlout_wgcrsoverview(const char *name, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s,  n(x, y)', 'Color', 'w');\n", name);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "lfs = 12;\n");
fprintf(dat, "alw = 1.5;\n");
fprintf(dat, "box on; hold on;\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(yend-ybeg)/(xend-xbeg) 0.95 0.95]);\n");
fprintf(dat, "xlabel(['y [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', lfs);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Draw refractive index patches\n");
fprintf(dat, "maxn = max(npn);\n");
fprintf(dat, "minn = min(npn);\n");
fprintf(dat, "for i=1:numnp\n");
fprintf(dat, "   pc = 1.0-(npn(i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "   patch([npy0(i) npy0(i) npy1(i) npy1(i)], ...\n");
fprintf(dat, "         [npx1(i) npx0(i) npx0(i) npx1(i)], ...\n");
fprintf(dat, "         [pc pc pc], 'LineStyle', 'none');\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Permittivity boundaries\n");
fprintf(dat, "lw = 1.5; ls = 'k-';\n");
	mlout_geooutline(dat);
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([ybeg ybeg yend yend ybeg], ...\n");
fprintf(dat, "     [xend xbeg xbeg xend xend], 'k-', 'LineWidth', 1.5);\n");
fprintf(dat, "\n");
	return;
}

/* ------------------------------------------------------------------- 
   routines related to the VEIMS solver, channel waveguides 
   ------------------------------------------------------------------- */

/* waveguide geometry and field extrema */
void mlout_geo(FILE *dat, WgCrs wg, double minf, double maxf)
{
	int i, j;
	fprintf(dat, "%% Waveguide geometry\n");
	fprintf(dat, "nlY = %d;\n", wg.ny);
	fprintf(dat, "bdY = [\n");
	for(j=0; j<=wg.ny; ++j) fprintf(dat, "%g ...\n", wg.hy(j));
	fprintf(dat, "];\n");
	fprintf(dat, "nlX = %d;\n", wg.nx);
	fprintf(dat, "bdX = [\n");
	for(i=0; i<=wg.nx; ++i) fprintf(dat, "%g ...\n", wg.hx(i));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Refractive indices\n");
	fprintf(dat, "n = [\n");
	for(i=0; i<=wg.nx+1; ++i) 
	{
		for(j=0; j<=wg.ny; ++j) 
			fprintf(dat, "%g ...\n", wg.n(i, j));
		fprintf(dat, "%g ;\n", wg.n(i, wg.ny+1));
	}
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Bounds on field values\n");
	fprintf(dat, "minf = %g;\n", minf);
	fprintf(dat, "maxf = %g;\n", maxf);
	fprintf(dat, "\n");
	return;
}

/* waveguide geometry and field extrema for interference visualization */
void mlout_geoyz(FILE *dat, WgCrs wg, double minf, double maxf)
{
	int i, j;
	double r;
	fprintf(dat, "%% Waveguide geometry\n");
	fprintf(dat, "nlY = %d;\n", wg.ny);
	fprintf(dat, "bdY = [\n");
	for(j=0; j<=wg.ny; ++j) fprintf(dat, "%g ...\n", wg.hy(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Refractive indices\n");
	fprintf(dat, "n = [\n");
	for(j=0; j<=wg.ny+1; ++j) 
	{
		r = 0.0;
		for(i=0; i<=wg.nx+1; ++i) r += wg.n(i, j);
		r = r/(wg.nx+2);
		fprintf(dat, "%g ...\n", r);
	}
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Bounds on field values\n");
	fprintf(dat, "minf = %g;\n", minf);
	fprintf(dat, "maxf = %g;\n", maxf);
	fprintf(dat, "\n");
	return;
}

/* three segment coupler geometry and field extrema for interference 
   visualization, input and output waveguides included */
#define BDDIST 1.0e-8
void mlout_iopropgeo(FILE *dat, WgCrs iwg, WgCrs pwg, WgCrs owg,
                     double l, double minf, double maxf)
{
	int i, j, k;
	double r;
	int p, pi, pp, po;
	double y, iy, py, oy;
	int xi;
	fprintf(dat, "%% Three segment coupler geometry\n");
	fprintf(dat, "nlZ = 1;\n");
	fprintf(dat, "bdZ = [ 0.0, %g];\n", l);
	Dvector bdy(iwg.ny+pwg.ny+owg.ny+3);
	p = pi = pp = po = 0; 
	while(pi <= iwg.ny || pp <= pwg.ny || po <= owg.ny)
	{
		if(pi <= iwg.ny) iy = iwg.hy(pi);
		else             iy = XYplane.y1;
		if(pp <= pwg.ny) py = pwg.hy(pp);
		else             py = XYplane.y1;
		if(po <= owg.ny) oy = owg.hy(po);
		else             oy = XYplane.y1;
		y = iy; 
		if(py < y) y = py;
		if(oy < y) y = oy;
		if(fabs(y-iy) < BDDIST) ++pi;
		if(fabs(y-py) < BDDIST) ++pp;
		if(fabs(y-oy) < BDDIST) ++po;
		bdy(p) = y;		
		++p;
	}
	--p;
	fprintf(dat, "nlY = %d;\n", p);
	fprintf(dat, "bdY = [\n");
	for(j=0; j<=p; ++j) fprintf(dat, "%g ...\n", bdy(j));
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Refractive indices\n");
	fprintf(dat, "n = [\n");
	for(k=0; k<=p+1; ++k) 
	{
		if(k==0) y = bdy(0)-1.0;		
		else
		{
			if(k==p+1) y = bdy(p)+1.0;
			else y = 0.5*(bdy(k)+bdy(k-1));
		}
		r = 0.0; iwg.rectidx(0.0, y, xi, j);
		for(i=0; i<=iwg.nx+1; ++i) r += iwg.n(i, j);
		r = r/(iwg.nx+2);
		fprintf(dat, "%g ...\n", r);
		r = 0.0; pwg.rectidx(0.0, y, xi, j);
		for(i=0; i<=pwg.nx+1; ++i) r += pwg.n(i, j);
		r = r/(pwg.nx+2);
		fprintf(dat, "%g ...\n", r);
		r = 0.0; owg.rectidx(0.0, y, xi, j);
		for(i=0; i<=owg.nx+1; ++i) r += owg.n(i, j);
		r = r/(pwg.nx+2);
		fprintf(dat, "%g ;\n", r);
	}
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	fprintf(dat, "%% Bounds on field values\n");
	fprintf(dat, "minf = %g;\n", minf);
	fprintf(dat, "maxf = %g;\n", maxf);
	fprintf(dat, "\n");
	return;
}

/* output mesh, x-y-plane */
void mlout_meshxy(FILE *dat, Rect disp, int npx, int npy)
{
	fprintf(dat, "%% Mesh information\n");
	fprintf(dat, "ybeg = %g;\n", disp.y0);
	fprintf(dat, "yend = %g;\n", disp.y1);
	fprintf(dat, "numy = %d;\n", npy);
	if(disp.y0==disp.y1)
		fprintf(dat, "y = %g*ones(numy, 1);\n", disp.y0);
	else
		fprintf(dat, "y = linspace(ybeg, yend, numy);\n");
	
	fprintf(dat, "xbeg = %g;\n", disp.x0);
	fprintf(dat, "xend = %g;\n", disp.x1);
	fprintf(dat, "numx = %d;\n", npx);
	if(disp.x0==disp.x1)
		fprintf(dat, "x = %g*ones(numx, 1);\n", disp.x0);
	else
		fprintf(dat, "x = linspace(xbeg, xend, numx);\n");
	fprintf(dat, "\n");
	return;
}

/* output mesh, y-z-plane */
void mlout_meshyz(FILE *dat, double ybeg, double yend, int npy,
                             double zbeg, double zend, int npz)
{
	fprintf(dat, "%% Mesh information\n");
	fprintf(dat, "ybeg = %g;\n", ybeg);
	fprintf(dat, "yend = %g;\n", yend);
	fprintf(dat, "numy = %d;\n", npy);
	if(ybeg==yend)
		fprintf(dat, "y = %g*ones(numy, 1);\n", ybeg);
	else
		fprintf(dat, "y = linspace(ybeg, yend, numy);\n");
	
	fprintf(dat, "zbeg = %g;\n", zbeg);
	fprintf(dat, "zend = %g;\n", zend);
	fprintf(dat, "numz = %d;\n", npz);
	if(zbeg==zend)
		fprintf(dat, "z = %g*ones(numz, 1);\n", zbeg);
	else
		fprintf(dat, "z = linspace(zbeg, zend, numz);\n");
	fprintf(dat, "\n");
	return;
}

/* field section */
void mlout_sec(FILE *dat, double x0, double y0, double x1, double y1, 
                     int np, Fcomp cp, char c2, char c3, Dvector fld)
{
	double dx, dy;
	int j;

	dx=(x1-x0)/(np-1); 
	dy=(y1-y0)/(np-1);
	fprintf(dat, "%c%c%c%c = [\n", 
                     fldchr(cp), cpchr(cp), c2, c3);
	for(j=0; j<=np-1; ++j) fprintf(dat, "%g\n", fld(j));
	fprintf(dat, "];\n");
	fprintf(dat, "xp%c%c = [\n", c2, c3);
	for(j=0; j<=np-1; ++j) fprintf(dat, "%g\n", x0+j*dx); 
	fprintf(dat, "];\n");
	fprintf(dat, "yp%c%c = [\n", c2, c3);
	for(j=0; j<=np-1; ++j) fprintf(dat, "%g\n", y0+j*dy);
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	return;
}

/* field section, 1D*/
void mlout_sec1D(FILE *dat,  
		 Fcomp cp, char c2, char c3, char c4, char c5,
                 int np, Dvector fld, Dvector pos)
{
	int j;

	fprintf(dat, "%c%c%c%c%c%c = [\n", 
                     fldchr(cp), cpchr(cp), c2, c3, c4, c5);
	for(j=0; j<=np-1; ++j) fprintf(dat, "%g\n", fld(j));
	fprintf(dat, "];\n");
	fprintf(dat, "pos%c%c%c%c = [\n", c2, c3, c4, c5);
	for(j=0; j<=np-1; ++j) fprintf(dat, "%g\n", pos(j)); 
	fprintf(dat, "];\n");
	fprintf(dat, "\n");
	return;
}

/* contour plot */
void mlout_contour(char *name, FILE *dat, Fcomp cp, Afo foa,
                   const char *vlab, const char *hlab)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "hold on;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "axis equal;\n");
fprintf(dat, "alw = 1.0;\n");
fprintf(dat, "bfs = abs((xend-xbeg)/(yend-ybeg));\n");
fprintf(dat, "if bfs>1 bfs=1/bfs; end;\n");
fprintf(dat, "bfs=11+5*bfs;\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', bfs, ...\n");
fprintf(dat, "	 'LineWidth', alw, ...\n");
fprintf(dat, "	 'TickDir', 'out');\n");
fprintf(dat, "xlabel(['%s [' char(956) 'm]'], 'FontSize', bfs);\n", hlab);
fprintf(dat, "ylabel(['%s [' char(956) 'm]'], 'FontSize', bfs);\n", vlab);
	switch(foa)
	{
		case MOD:  
fprintf(dat, "title('|%c_%c|', 'FontSize', bfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case ORG: 
fprintf(dat, "title('%c_%c', 'FontSize', bfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case SQR:  
fprintf(dat, "title('|%c_%c|^2', 'FontSize', bfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case REP:  
fprintf(dat, "title('Re %c_%c', 'FontSize', bfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case IMP:  
fprintf(dat, "title('Im %c_%c', 'FontSize', bfs);\n",
             fldchr(cp), cpchr(cp));
			break;
	}
fprintf(dat, "\n");
fprintf(dat, "%% Background: shading according to waveguide geometry\n");
fprintf(dat, "maxn = max(max(n));\n");
fprintf(dat, "minn = min(min(n));\n");
fprintf(dat, "for j=1:nlX+2\n");
fprintf(dat, "   if j==1 x0=xbeg;\n");
fprintf(dat, "   else    x0=bdX(j-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if j==nlX+2 x1=xend;\n");
fprintf(dat, "   else        x1=bdX(j);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for i=1:nlY+2\n");
fprintf(dat, "      if i==1 y0=ybeg;\n");
fprintf(dat, "      else    y0=bdY(i-1);\n");
fprintf(dat, "      end\n");
fprintf(dat, "      if i==nlY+2 y1=yend;\n");
fprintf(dat, "      else        y1=bdY(i);\n");
fprintf(dat, "      end\n");
fprintf(dat, "      pc = 1.0-(n(j, i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "      patch([y0 y0 y1 y1], [x1 x0 x0 x1], [pc pc pc], ...\n");
fprintf(dat, "        'LineStyle', 'none');\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add field contours\n");
	if(Mlop_Colour==YES) fprintf(dat, "cl = 'b-';\n");
	else                 fprintf(dat, "cl = 'k-';\n");
fprintf(dat, "[c, h] =contour(y, x, %c%c, 19, cl); \n", 
              fldchr(cp), cpchr(cp));
fprintf(dat, "set(h, 'LineWidth', 0.5);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Indicate permittivity boundaries\n");
fprintf(dat, "lw = 0.5;\n");
fprintf(dat, "for j=1:nlX+2\n");
fprintf(dat, "   if j==1 x0=xbeg;\n");
fprintf(dat, "   else    x0=bdX(j-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if j==nlX+2 x1=xend;\n");
fprintf(dat, "   else        x1=bdX(j);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for i=1:nlY+1\n");
fprintf(dat, "      yp=bdY(i); \n");
fprintf(dat, "      if n(j,i)~=n(j,i+1)\n");
fprintf(dat, "         plot([yp, yp], [x0, x1], 'k-', 'LineWidth', lw);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for i=1:nlY+2\n");
fprintf(dat, "   if i==1 y0=ybeg;\n");
fprintf(dat, "   else    y0=bdY(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nlY+2 y1=yend;\n");
fprintf(dat, "   else        y1=bdY(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for j=1:nlX+1\n");
fprintf(dat, "      xp=bdX(j); \n");
fprintf(dat, "      if n(j,i)~=n(j+1,i)\n");
fprintf(dat, "         plot([y0, y1], [xp, xp], 'k-', 'LineWidth', lw);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([ybeg ybeg yend yend ybeg], ...\n");
fprintf(dat, "     [xend xbeg xbeg xend xend], 'k-', 'LineWidth', alw);\n");
fprintf(dat, "\n");
	return;
}
void mlout_contour(char *name, FILE *dat, Fcomp cp, Afo foa)
{
	mlout_contour(name, dat, cp, foa, "x", "y");
	return;
}

/* surface plot */
void mlout_surface(char *name, FILE *dat, Fcomp cp, Afo foa,
                   const char *vlab, const char *hlab)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "\n");
fprintf(dat, "%% Adjust field aspect ratio to cross section extension\n");
fprintf(dat, "fld = %c%c;\n", fldchr(cp), cpchr(cp));
fprintf(dat, "xyzfac = min([yend-ybeg, xend-xbeg])*0.8;\n");
fprintf(dat, "%% uncomment to adjust plot to actual component\n");
fprintf(dat, "%% maxf = max(max(fld));\n");
fprintf(dat, "%% minf = min(min(fld));\n");
fprintf(dat, "\n");
fprintf(dat, "df   = maxf-minf;\n");
fprintf(dat, "fld  =  fld/df*xyzfac;\n");
fprintf(dat, "minf = minf/df*xyzfac;\n");
fprintf(dat, "maxf = maxf/df*xyzfac;\n");
fprintf(dat, "\n");
	if(Mlop_Colour==YES)
	{
fprintf(dat, "%% Choose colormap\n");
fprintf(dat, "black = zeros(64, 3);\n");
fprintf(dat, "blue = black;\n");
fprintf(dat, "blue(:,3) = 1;\n");
fprintf(dat, "colormap(blue);\n");
	}
	else
	{
fprintf(dat, "black = zeros(64, 3);\n");
fprintf(dat, "colormap(black);\n");
	}
fprintf(dat, "\n");
fprintf(dat, "%% Surface plot \n");
fprintf(dat, "mesh(y, x, fld, 'LineWidth', 0.1);\n");
fprintf(dat, "%% surf(y, x, fld, 'LineWidth', 0.1);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minf, maxf]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "set(gca, 'FontSize', 15);\n");
fprintf(dat, "xlabel(['%s [' char(956) 'm]'], 'FontSize', 17);\n", hlab);
fprintf(dat, "ylabel(['%s [' char(956) 'm]'], 'FontSize', 17);\n", vlab);
	switch(foa)
	{
		case MOD:  
fprintf(dat, "zlabel('|%c_%c| [a.u.]', 'FontSize', 17);\n",
             fldchr(cp), cpchr(cp));
			break;
		case ORG: 
fprintf(dat, "zlabel('%c_%c [a.u.]', 'FontSize', 17);\n",
             fldchr(cp), cpchr(cp));
			break;
		case SQR:  
fprintf(dat, "zlabel('|%c_%c|^2 [a.u.]', 'FontSize', 17);\n",
             fldchr(cp), cpchr(cp));
			break;
		case REP:  
fprintf(dat, "zlabel('Re %c_%c [a.u.]', 'FontSize', 17);\n",
             fldchr(cp), cpchr(cp));
			break;
		case IMP:  
fprintf(dat, "zlabel('Im %c_%c [a.u.]', 'FontSize', 17);\n",
             fldchr(cp), cpchr(cp));
			break;
	}
	switch(foa)
	{
		case MOD:  
fprintf(dat, "%% title('|%c_%c|', 'FontSize', 20);\n",
             fldchr(cp), cpchr(cp));
			break;
		case ORG: 
fprintf(dat, "%% title('%c_%c', 'FontSize', 20);\n",
             fldchr(cp), cpchr(cp));
			break;
		case SQR:  
fprintf(dat, "%% title('|%c_%c|^2', 'FontSize', 20);\n",
             fldchr(cp), cpchr(cp));
			break;
		case REP:  
fprintf(dat, "%% title('Re %c_%c', 'FontSize', 20);\n",
             fldchr(cp), cpchr(cp));
			break;
		case IMP:  
fprintf(dat, "%% title('Im %c_%c', 'FontSize', 20);\n",
             fldchr(cp), cpchr(cp));
			break;
	}
fprintf(dat, "%% axis off;\n");
fprintf(dat, "%% grid off;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "view(-130,20);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries \n");
fprintf(dat, "hold on;\n");
fprintf(dat, "lw = 2.5;\n");
fprintf(dat, "zp = maxf;\n");
fprintf(dat, "%% zp = minf;\n");
fprintf(dat, "lc = 'k';\n");
fprintf(dat, "for j=1:nlX+2\n");
fprintf(dat, "   if j==1 x0=xbeg;\n");
fprintf(dat, "   else    x0=bdX(j-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if j==nlX+2 x1=xend;\n");
fprintf(dat, "   else        x1=bdX(j);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for i=1:nlY+1\n");
fprintf(dat, "      yp=bdY(i); \n");
fprintf(dat, "      if n(j,i)~=n(j,i+1)\n");
fprintf(dat, "         line([yp, yp], [x0, x1], [zp, zp], ... \n");
fprintf(dat, "              'LineWidth', lw, 'Color', lc);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end   \n");
fprintf(dat, "for i=1:nlY+2\n");
fprintf(dat, "   if i==1 y0=ybeg;\n");
fprintf(dat, "   else    y0=bdY(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nlY+2 y1=yend;\n");
fprintf(dat, "   else        y1=bdY(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for j=1:nlX+1\n");
fprintf(dat, "      xp=bdX(j); \n");
fprintf(dat, "      if n(j,i)~=n(j+1,i)\n");
fprintf(dat, "   	   line([y0, y1], [xp, xp], [zp, zp], ... \n");
fprintf(dat, "                'LineWidth', lw, 'Color', lc);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}
void mlout_surface(char *name, FILE *dat, Fcomp cp, Afo foa)
{
	mlout_surface(name, dat, cp, foa, "x", "y");
	return;
}

/* intensity image */
void mlout_image(char *name, FILE *dat, Fcomp cp, Afo foa, 
                 double minf, double maxf,
                 const char *vlab, const char *hlab)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "minf = %.10g;\n", minf);
fprintf(dat, "maxf = %.10g;\n", maxf);
fprintf(dat, "%% minf = min(min(%c%c));\n", fldchr(cp), cpchr(cp));
fprintf(dat, "%% maxf = max(max(%c%c));\n", fldchr(cp), cpchr(cp));
fprintf(dat, "if minf >= maxf\n");
fprintf(dat, "   minf = 0.5*(maxf+minf)-0.5;\n");
fprintf(dat, "   maxf = minf+1.0;\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Intensity image\n");
fprintf(dat, "imagesc(y, x, %c%c, [minf, maxf]);\n", fldchr(cp), cpchr(cp));
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "lfs = 12;\n");
fprintf(dat, "alw = 1.5;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', lfs, ...\n");
fprintf(dat, "         'LineWidth', alw, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [(yend-ybeg)/(xend-xbeg) 0.95 0.95]);\n");
fprintf(dat, "xlabel(['%s [' char(956) 'm]'], 'FontSize', lfs);\n", hlab);
fprintf(dat, "ylabel(['%s [' char(956) 'm]'], 'FontSize', lfs);\n", vlab);
fprintf(dat, "\n");
	switch(foa)
	{
		case MOD:  
fprintf(dat, "title('|%c_%c|', 'FontSize', lfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case ORG: 
fprintf(dat, "title('%c_%c', 'FontSize', lfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case SQR:  
fprintf(dat, "title('|%c_%c|^2', 'FontSize', lfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case REP:  
fprintf(dat, "title('Re %c_%c', 'FontSize', lfs);\n",
             fldchr(cp), cpchr(cp));
			break;
		case IMP:  
fprintf(dat, "title('Im %c_%c', 'FontSize', lfs);\n",
             fldchr(cp), cpchr(cp));
			break;
	}
fprintf(dat, "\n");
	mlout_sketchinit(dat);
	mlout_geooutline(dat);
	return;
}
void mlout_image(char *name, FILE *dat, Fcomp cp, Afo foa, 
                 double minf, double maxf)
{
	mlout_image(name, dat, cp, foa, minf, maxf, "x", "y");
	return;
}

/* six component surface plot */  
void mlout_acmfile(char *name, FILE *dat, 
		   double minE, double maxE, double minH, double maxH)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "\n");
fprintf(dat, "%% Field boundary values\n");
fprintf(dat, "minE = %g;\n", minE);
fprintf(dat, "maxE = %g;\n", maxE);
fprintf(dat, "minH = %g;\n", minH);
fprintf(dat, "maxH = %g;\n", maxH);
fprintf(dat, "\n");
fprintf(dat, "%% Adjust field aspect ratio to cross section extension\n");
fprintf(dat, "s = 1; %% set to -1 to invert plot\n");
fprintf(dat, "xyzfac = min([yend-ybeg, xend-xbeg])*0.8;\n");
fprintf(dat, "dE = maxE-minE;\n");
fprintf(dat, "Ex = s*Ex/dE*xyzfac;\n");
fprintf(dat, "Ey = s*Ey/dE*xyzfac;\n");
fprintf(dat, "Ez = s*Ez/dE*xyzfac;\n");
fprintf(dat, "minE = s*minE/dE*xyzfac;\n");
fprintf(dat, "maxE = s*maxE/dE*xyzfac;\n");
fprintf(dat, "if minE > maxE tE=minE; minE=maxE; maxE=tE; end\n");
fprintf(dat, "dH = maxH-minH;\n");
fprintf(dat, "Hx = s*Hx/dH*xyzfac;\n");
fprintf(dat, "Hy = s*Hy/dH*xyzfac;\n");
fprintf(dat, "Hz = s*Hz/dH*xyzfac;\n");
fprintf(dat, "minH = s*minH/dH*xyzfac;\n");
fprintf(dat, "maxH = s*maxH/dH*xyzfac;\n");
fprintf(dat, "if minH > maxH tH=minH; minH=maxH; maxH=tH; end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Title size\n");
fprintf(dat, "ts = 15;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Line width\n");
fprintf(dat, "lw = 0.1;\n");
fprintf(dat, "\n");
	if(Mlop_Colour==YES)
	{
fprintf(dat, "%% Choose colormap\n");
fprintf(dat, "black = zeros(64, 3);\n");
fprintf(dat, "blue = black;\n");
fprintf(dat, "blue(:,3) = 1;\n");
fprintf(dat, "colormap(blue);\n");
	}
	else
	{
fprintf(dat, "black = zeros(64, 3);\n");
fprintf(dat, "colormap(black);\n");
	}
fprintf(dat, "\n");
fprintf(dat, "%% E_x: top row, left\n");
fprintf(dat, "subplot('position', [0.01 0.51 0.31 0.49])\n");
fprintf(dat, "mesh(y, x, Ex, 'LineWidth', lw);\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minE, maxE]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "title('E_x', 'FontSize', ts);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "\n");
fprintf(dat, "%% E_y: top row, center\n");
fprintf(dat, "subplot('position', [0.34 0.51 0.31 0.49])\n");
fprintf(dat, "mesh(y, x, Ey, 'LineWidth', lw);\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minE, maxE]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "title('E_y', 'FontSize', ts);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "\n");
fprintf(dat, "%% E_z: top row, right\n");
fprintf(dat, "subplot('position', [0.67 0.51 0.31 0.49])\n");
fprintf(dat, "mesh(y, x, Ez, 'LineWidth', lw);\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minE, maxE]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "title('E_z', 'FontSize', ts);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "\n");
fprintf(dat, "%% H_x: bottom row, left\n");
fprintf(dat, "subplot('position', [0.01 0.01 0.31 0.49])\n");
fprintf(dat, "mesh(y, x, Hx, 'LineWidth', lw);\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minH, maxH]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "title('H_x', 'FontSize', ts);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "\n");
fprintf(dat, "%% H_y: bottom row, center\n");
fprintf(dat, "subplot('position', [0.34 0.01 0.31 0.49])\n");
fprintf(dat, "mesh(y, x, Hy, 'LineWidth', lw);\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minH, maxH]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "title('H_y', 'FontSize', ts);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "\n");
fprintf(dat, "%% H_z: bottom row, center\n");
fprintf(dat, "subplot('position', [0.67 0.01 0.31 0.49])\n");
fprintf(dat, "mesh(y, x, Hz, 'LineWidth', lw);\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minH, maxH]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1]);\n");
fprintf(dat, "title('H_z', 'FontSize', ts);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "\n");
	return;
}

/* mode section plot geometry  */
void mlout_Lsecgeo(FILE *dat, double x0, double y0, double x1, double y1,
                   int nbd, Dvector bd, Dvector ri)
{
	int i;
fprintf(dat, "%% section location\n");
fprintf(dat, "x0 = %g;\n", x0);
fprintf(dat, "x1 = %g;\n", x1);
fprintf(dat, "y0 = %g;\n", y0);
fprintf(dat, "y1 = %g;\n", y1);
fprintf(dat, "\n");
fprintf(dat, "%% Crossed boundaries \n");
fprintf(dat, "nbd = %d;\n", nbd);
fprintf(dat, "bd = [ \n");
	for(i=0; i<=nbd-1; ++i)
	{
fprintf(dat, "%g\n", bd(i));
	}
fprintf(dat, "];\n");
fprintf(dat, "ri = [ \n");
	for(i=0; i<=nbd; ++i)
	{
fprintf(dat, "%g\n", ri(i));
	}
fprintf(dat, "];\n");
fprintf(dat, "\n");
	return;
}

/* mode section plot, plot commands */
void mlout_Lsecplot(char *name, FILE *dat, char ori, double minf, double maxf,
                    Fcomp cp, int nbd, int nsec)
{
	int i;
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "hold on;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
fprintf(dat, "%% uncomment to adjust axes and labels to current section\n");
fprintf(dat, "%% minf = %g;\n", minf);
fprintf(dat, "%% maxf = %g;\n", maxf);
	if(ori=='h') 
	{
fprintf(dat, "   a0 = y0;\n");
fprintf(dat, "   a1 = y1; \n");
fprintf(dat, "   xlabel(['y [' char(956) 'm]'], 'FontSize', 17);\n");
	}
	else
	{
fprintf(dat, "   a0 = x0;\n");
fprintf(dat, "   a1 = x1; \n");
fprintf(dat, "   xlabel(['x [' char(956) 'm]'], 'FontSize', 17);\n");
	}
fprintf(dat, "axis([a0, a1, minf, maxf]);\n");
fprintf(dat, "alw = 1.0;\n");
fprintf(dat, "set(gca, 'FontSize', 15, ...\n");
fprintf(dat, "	 'LineWidth', alw, ...\n");
fprintf(dat, "	 'TickDir', 'out');\n");
fprintf(dat, "ylabel('%c_%c', 'FontSize', 17);\n", fldchr(cp), cpchr(cp));
fprintf(dat, "\n");
	if(nbd > 0)
	{
fprintf(dat, "%% Background: shading according to waveguide geometry\n");
fprintf(dat, "maxn = max(max(n));\n");
fprintf(dat, "minn = min(min(n));\n");
fprintf(dat, "for i=1:nbd+1\n");
fprintf(dat, "   if i==1 b0=a0;\n");
fprintf(dat, "   else    b0=bd(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nbd+1 b1=a1;\n");
fprintf(dat, "   else        b1=bd(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   pc = 1.0-(ri(i)-minn)/(maxn-minn)*0.2;\n");
fprintf(dat, "   patch([b0 b0 b1 b1], [maxf minf minf maxf], [pc pc pc], ...\n");
fprintf(dat, "           'LineStyle', 'none');\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	}
fprintf(dat, "%% Add field values\n");
fprintf(dat, "lw = 2.0;\n");
	if(Mlop_Colour==YES)
	{
fprintf(dat, "ls = 'b-';\n");
	}
	else
	{
fprintf(dat, "ls = 'k-';\n");
	}
	for(i=0; i<=nsec-1; ++i)
	{
fprintf(dat, "plot(pos%c%c, %c%c%c%c, ls, 'LineWidth', lw);\n", 
	     dig10(i), dig1(i), 
             fldchr(cp), cpchr(cp), dig10(i), dig1(i));
	}
fprintf(dat, "\n");
fprintf(dat, "%% Include further field section data\n");
fprintf(dat, "%% _______V\n");
fprintf(dat, "\n");
	if(nbd > 0)
	{
fprintf(dat, "%% Indicate permittivity boundaries\n");
fprintf(dat, "for i=1:nbd\n");
fprintf(dat, "   plot([bd(i), bd(i)], [minf maxf], 'k-', 'LineWidth', 0.5);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Renew frame\n");
fprintf(dat, "plot([a0 a0 a1 a1 a0], ...\n");
fprintf(dat, "     [maxf minf minf maxf maxf], 'k-', 'LineWidth', alw);\n");
fprintf(dat, "\n");
	}
	return;
}

/* fancy plot */
void mlout_fancy(char *name, FILE *dat, Fcomp cp, int numc)
{
	int j;
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "\n");
fprintf(dat, "%% Adjust field aspect ratio to cross section extension\n");
fprintf(dat, "s = 1; %% set to -1 to invert the plot \n");
fprintf(dat, "fld = %c%c;\n", fldchr(cp), cpchr(cp));
fprintf(dat, "xyzfac = min([yend-ybeg, xend-xbeg])*0.5;\n");
fprintf(dat, "%% uncomment to adjust plot to actual field\n");
fprintf(dat, "%% maxf = max(max(fld));\n");
fprintf(dat, "%% minf = min(min(fld));\n");
fprintf(dat, "df   = maxf-minf;\n");
fprintf(dat, "fld  =  s*fld/df*xyzfac;\n");
fprintf(dat, "minf = s*minf/df*xyzfac;\n");
fprintf(dat, "maxf = s*maxf/df*xyzfac;\n");
fprintf(dat, "if minf > maxf tf=minf; minf=maxf; maxf=tf; end\n");
fprintf(dat, "\n");
	mlout_magcolormap(dat);
fprintf(dat, "\n");
fprintf(dat, "%% Surface plot \n");
fprintf(dat, "sf = surf(y, x, fld);\n");
fprintf(dat, "set(sf, 'LineStyle', 'none');\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend, minf, maxf]);\n");
fprintf(dat, "set(gca, 'DataAspectRatio', [1, 1, 1], ...\n");
fprintf(dat, "         'Projection', 'Perspective', ...\n");
fprintf(dat, "         'Units', 'normalized', ...\n");
fprintf(dat, "         'Position', [0, 0, 1, 1]);\n");
fprintf(dat, "axis off;\n");
fprintf(dat, "grid off;\n");
fprintf(dat, "view(-130,20);\n");
fprintf(dat, "caxis([-maxf maxf]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Lighting\n");
fprintf(dat, "lighting phong;\n");
fprintf(dat, "light('Position', [(yend+ybeg)/2, xbeg-2*(xend-xbeg), maxf+2*(maxf-minf)]);\n");
fprintf(dat, "light('Position', [(yend+ybeg)/2, xend+2*(xend-xbeg), maxf+2*(maxf-minf)]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries \n");
fprintf(dat, "hold on;\n");
	if(Mlop_Print==YES) fprintf(dat, "lw = 3.0;\n");
	else                fprintf(dat, "lw = 3.0;\n");
fprintf(dat, "lc = 'k';\n");
fprintf(dat, "lift = s*df/200.0;\n");
	for(j=0; j<=numc-1; ++j)
	{
fprintf(dat, 
"line(yp%c%c, xp%c%c, s*(%c%c%c%c+lift)/df*xyzfac, ...\n",
             dig10(j), dig1(j), dig10(j), dig1(j), 
             fldchr(cp), cpchr(cp), dig10(j), dig1(j));
fprintf(dat, "   'LineWidth', lw, 'Color', lc);\n");
	}
fprintf(dat, "\n");
	return;
}

/* for interference data: show the z position */
void mlout_annotatezpos(FILE *dat, double rz)
{
fprintf(dat, "%% show z-position\n");
fprintf(dat, "text(yend-(yend-ybeg)/40, xbeg+(xend-xbeg)/20, ...\n");
fprintf(dat, "   ['z=%g' char(956) 'm]'], ...\n", rz);
fprintf(dat, "   'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "   'VerticalAlignment', 'bottom', 'FontSize', bfs);\n");
	return;
}

/* .m file title comment */
void mlout_title(FILE *dat, char *name, char *intro)
{
fprintf(dat, "%% %s\n", name);
fprintf(dat, "%% %s\n", intro);
fprintf(dat, "\n");
return;
}
 
/* interference field animation */
void mlout_play(FILE *dat, char *name, int numz)
{
	int i;
fprintf(dat, "%% first frame\n");
	name[6] = '0';
	name[7] = '0';
fprintf(dat, "%s;\n", name);
fprintf(dat, "%% reserve memory\n");
fprintf(dat, "M = moviein(%d);\n", numz);
fprintf(dat, "%% store frames\n");
fprintf(dat, "M(:,1) = getframe;\n");
	for(i=1; i<=numz-1; ++i)
	{
fprintf(dat, "close;\n");
		name[6] = dig10(i);
		name[7] = dig1(i);
fprintf(dat, "%s;\n", name);
fprintf(dat, "M(:,%d) = getframe;\n", i+1);
	}
fprintf(dat, "set(gcf, 'Name', 'Mode Interference');\n");
fprintf(dat, "\n");
fprintf(dat, "%% ... show it. \n");
fprintf(dat, "movie(M, 10)\n");
	return;
}

/* interference intensity image */
void mlout_propimage(char *name, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
	if(Mlop_Colour==YES)
	{
fprintf(dat, "%% Choose colormap\n");
fprintf(dat, "blue = zeros(64, 3);\n");
fprintf(dat, "blue(:,1) = linspace(0,1,64)';\n");
fprintf(dat, "blue(:,2) = blue(:,1);\n");
fprintf(dat, "blue(:,3) = 1;\n");
fprintf(dat, "colormap(blue);\n");
	}
	else
	{
fprintf(dat, "colormap(gray);\n");
	}
fprintf(dat, "\n");
fprintf(dat, "%% adjust intensity range to actual field (comment for a composite plot)\n");
fprintf(dat, "maxf = max(max(Sz));\n");
fprintf(dat, "minf = min(min(Sz));\n");
fprintf(dat, "\n");
fprintf(dat, "%% Intensity image\n");
fprintf(dat, "imagesc(z, y, Sz, [minf, maxf]);\n"); 
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "box on;\n");
fprintf(dat, "axis([zbeg, zend, ybeg, yend]);\n");
fprintf(dat, "set(gca, 'FontSize', 10, ...\n");
fprintf(dat, "         'LineWidth', 1.5, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [3.5 1 1]);\n");
fprintf(dat, "xlabel(['z [' char(956) 'm]'], 'FontSize', 12);\n");
fprintf(dat, "ylabel(['y [' char(956) 'm]'], 'FontSize', 12);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "lw = 1.0;\n");
	if(Mlop_Colour==YES)
fprintf(dat, "ls = 'k-';\n");
	else
fprintf(dat, "ls = 'w-';\n");
fprintf(dat, "for i=1:nlY+1\n");
fprintf(dat, "   yp=bdY(i); \n");
fprintf(dat, "   if n(i)~=n(i+1)\n");
fprintf(dat, "      plot([zbeg, zend], [yp, yp], ls, 'LineWidth', lw);\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* three segment coupler, interference intensity image */
void mlout_iopropimage(char *name, FILE *dat)
{
fprintf(dat, "%% Initialize the plot\n");
fprintf(dat, "figure('NumberTitle', 'off', 'Name', '%s', 'Color', 'w');\n", name);
fprintf(dat, "view(2);\n");
fprintf(dat, "\n");
	if(Mlop_Colour==YES)
	{
fprintf(dat, "%% Choose colormap\n");
fprintf(dat, "blue = zeros(64, 3);\n");
fprintf(dat, "blue(:,1) = linspace(0,1,64)';\n");
fprintf(dat, "blue(:,2) = blue(:,1);\n");
fprintf(dat, "blue(:,3) = 1;\n");
fprintf(dat, "colormap(blue);\n");
	}
	else
	{
fprintf(dat, "colormap(gray);\n");
	}
fprintf(dat, "\n");
fprintf(dat, "%% adjust intensity range to actual field (comment for a composite plot)\n");
fprintf(dat, "maxf = max(max(Sz));\n");
fprintf(dat, "minf = min(min(Sz));\n");
fprintf(dat, "\n");
fprintf(dat, "%% Intensity image\n");
fprintf(dat, "imagesc(z, y, Sz, [minf, maxf]);\n"); 
fprintf(dat, "\n");
fprintf(dat, "%% Adjust axes and labels\n");
fprintf(dat, "box on;\n");
fprintf(dat, "axis([zbeg, zend, ybeg, yend]);\n");
fprintf(dat, "set(gca, 'FontSize', 10, ...\n");
fprintf(dat, "         'LineWidth', 1.5, ...\n");
fprintf(dat, "         'Clipping', 'off', ...\n");
fprintf(dat, "         'YDir', 'normal', ...\n");
fprintf(dat, "         'TickDir', 'out', ...\n");
fprintf(dat, "         'PlotBoxAspectRatio', [3.5 1 1]);\n");
fprintf(dat, "xlabel(['z [' char(956) 'm]'], 'FontSize', 12);\n");
fprintf(dat, "ylabel(['y [' char(956) 'm]'], 'FontSize', 12);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "lw = 1.0;\n");
	if(Mlop_Colour==YES)
fprintf(dat, "ls = 'k-';\n");
	else
fprintf(dat, "ls = 'w-';\n");
fprintf(dat, "for j=1:nlY+2\n");
fprintf(dat, "   if j==1 y0=ybeg;\n");
fprintf(dat, "   else    y0=bdY(j-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if j==nlY+2 y1=yend;\n");
fprintf(dat, "   else        y1=bdY(j);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for i=1:nlZ+1\n");
fprintf(dat, "      zp=bdZ(i); \n");
fprintf(dat, "      if n(j,i)~=n(j,i+1)\n");
fprintf(dat, "         plot([zp, zp], [y0, y1], ls, 'LineWidth', lw);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for i=1:nlZ+2\n");
fprintf(dat, "   if i==1 z0=zbeg;\n");
fprintf(dat, "   else    z0=bdZ(i-1);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   if i==nlZ+2 z1=zend;\n");
fprintf(dat, "   else        z1=bdZ(i);\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for j=1:nlY+1\n");
fprintf(dat, "      yp=bdY(j); \n");
fprintf(dat, "      if n(j,i)~=n(j+1,i)\n");
fprintf(dat, "         plot([z0, z1], [yp, yp], ls, 'LineWidth', lw);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* print commands */
void mlout_print(FILE *dat, char *name, char dev)
{
	if(Mlop_Print != YES) 
	{
		switch(dev)
		{
		case 'p':
fprintf(dat, "%%%% Print figure to file\n");
fprintf(dat, "%% set(gcf, 'InvertHardcopy', 'off');\n");
fprintf(dat, "%% print -dpng -r200 %s \n", name);
			break;
		default:
fprintf(dat, "%%%% Write encapsulated postscript file\n");
fprintf(dat, "%% set(gcf, 'InvertHardcopy', 'off');\n");
fprintf(dat, "%% print -depsc2 -r600 %s\n", name);
fprintf(dat, "\n");
			break;
		}
		return;	
	}
	switch(dev)
	{
	case 'p':
fprintf(dat, "%% Print figure to file\n");
fprintf(dat, "set(gcf, 'InvertHardcopy', 'off');\n");
fprintf(dat, "print -dpng -r200 %s \n", name);
		break;
	default:
fprintf(dat, "%% Write encapsulated postscript file\n");
fprintf(dat, "set(gcf, 'InvertHardcopy', 'off');\n");
fprintf(dat, "print -depsc2 -r600 %s\n", name);
fprintf(dat, "\n");
		break;
	}
	return;	
}
/* ------------------------------------------------------------------------ */ 

/* field viewer .m-file, 3D mode profiles: top, title and globals */
void mlout_viewertopxy(FILE *dat, const char *name, 
                       Polarization p, double lambda)
{
fprintf(dat, "%% %s.m\n", name);
fprintf(dat, "%% Metric VEIMS mode profile, viewer GUI\n");
fprintf(dat, "function %s\n", name);
fprintf(dat, "\n");
fprintf(dat, "close all; clear global;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Polarisation\n");
fprintf(dat, "pol = 'T%c';\n", polCHR(p));
fprintf(dat, "%% Wavelength\n");
fprintf(dat, "lambda = %.10g;\n", lambda);
fprintf(dat, "\n");
fprintf(dat, "%% Field and structure data\n");
fprintf(dat, "global nvd vdbeg vdend vdpos nhd hdbeg hdend hdpos n;\n");
fprintf(dat, "global xbeg xend numx x ybeg yend numy y;\n");
fprintf(dat, "global Ex_re Ex_im Ey_re Ey_im Ez_re Ez_im;\n");
fprintf(dat, "global Hx_re Hx_im Hy_re Hy_im Hz_re Hz_im;\n");
fprintf(dat, "\n");
return;
}
void mlout_viewertopxy(FILE *dat, const char *name, double lambda)
{
fprintf(dat, "%% %s.m\n", name);
fprintf(dat, "%% Metric VEIMS mode profile, viewer GUI\n");
fprintf(dat, "function %s\n", name);
fprintf(dat, "\n");
fprintf(dat, "close all; clear global;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Polarisation\n");
fprintf(dat, "pol = 'VEC';\n");
fprintf(dat, "%% Wavelength\n");
fprintf(dat, "lambda = %.10g;\n", lambda);
fprintf(dat, "\n");
fprintf(dat, "%% Field and structure data\n");
fprintf(dat, "global nvd vdbeg vdend vdpos nhd hdbeg hdend hdpos n;\n");
fprintf(dat, "global xbeg xend numx x ybeg yend numy y;\n");
fprintf(dat, "global Ex_re Ex_im Ey_re Ey_im Ez_re Ez_im;\n");
fprintf(dat, "global Hx_re Hx_im Hy_re Hy_im Hz_re Hz_im;\n");
fprintf(dat, "\n");
return;
}

/* mode profile viewer .m-file: viewer commands */
void mlout_fldviewerxy(FILE *dat, const char *name)
{
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% Viewer commands \n");
fprintf(dat, "\n");
fprintf(dat, "clear Ex; clear Ey; clear Ez; clear Hx; clear Hy; clear Hz;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Poynting vector\n");
fprintf(dat, "global Sz;\n");
fprintf(dat, "Sz=0.5*real((Ex_re+i*Ex_im).*(Hy_re-i*Hy_im)-(Ey_re+i*Ey_im).*(Hx_re-i*Hx_im)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Energy density of the optical field\n");
fprintf(dat, "global W;\n");
fprintf(dat, "mu0 = 1.2566370614E+3; ep0 = 8.854187817E-3;\n");
fprintf(dat, "epsr = sign(n).*n.*n;\n");
fprintf(dat, "W = 0.25*(mu0*(Hx_re.*Hx_re+Hx_im.*Hx_im+Hy_re.*Hy_re+Hy_im.*Hy_im ...\n");
fprintf(dat, "              +Hz_re.*Hz_re+Hz_im.*Hz_im ) ...\n");
fprintf(dat, "         +ep0.*epsr.*(Ex_re.*Ex_re+Ex_im.*Ex_im+Ey_re.*Ey_re+Ey_im.*Ey_im ...\n");
fprintf(dat, "                     +Ez_re.*Ez_re+Ez_im.*Ez_im ));\n");
fprintf(dat, "\n");
fprintf(dat, "%% Maximum field levels\n");
fprintf(dat, "global max_E max_H max_S max_W min_n max_n;\n");
fprintf(dat, "max_E = max(max(Ex_re.*Ex_re+Ex_im.*Ex_im));\n");
fprintf(dat, "m = max(max(Ey_re.*Ey_re+Ey_im.*Ey_im));\n");
fprintf(dat, "if m>max_E max_E = m; end\n");
fprintf(dat, "m = max(max(Ez_re.*Ez_re+Ez_im.*Ez_im));\n");
fprintf(dat, "if m>max_E max_E = m; end\n");
fprintf(dat, "max_E = sqrt(max_E);\n");
fprintf(dat, "max_H = max(max(Hx_re.*Hx_re+Hx_im.*Hx_im));\n");
fprintf(dat, "m = max(max(Hy_re.*Hy_re+Hy_im.*Hy_im));\n");
fprintf(dat, "if m>max_H max_H = m; end\n");
fprintf(dat, "m = max(max(Hz_re.*Hz_re+Hz_im.*Hz_im));\n");
fprintf(dat, "if m>max_H max_H = m; end\n");
fprintf(dat, "max_H = sqrt(max_H);\n");
fprintf(dat, "max_S = max(max(abs(Sz)));\n");
fprintf(dat, "max_W = max(max(abs(W)));\n");
fprintf(dat, "max_n = max(max(n));\n");
fprintf(dat, "min_n = min(min(n));\n");
fprintf(dat, "if min_n >= max_n max_n = min_n+1; end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Linear colormap\n");
fprintf(dat, "global lin_blue;\n");
fprintf(dat, "mi = 32; mlv = 0.2; lin_blue = zeros(64, 3);\n");
fprintf(dat, "lin_blue(1:mi, 3) = linspace(0, 1,   mi)';\n");
fprintf(dat, "lin_blue(1:mi, 2) = linspace(0, mlv, mi)';\n");
fprintf(dat, "lin_blue(1:mi, 1) = linspace(0, mlv, mi)';\n");
fprintf(dat, "lin_blue(mi+1:64, 3) = 1;\n");
fprintf(dat, "lin_blue(mi+1:64, 2) = linspace(mlv, 1, 64-mi)';\n");
fprintf(dat, "lin_blue(mi+1:64, 1) = linspace(mlv, 1, 64-mi)';\n");
fprintf(dat, "\n");
fprintf(dat, "%% Colormap: magnified zero level\n");
fprintf(dat, "global mag_blue;\n");
fprintf(dat, "mag_blue = zeros(256, 3);\n");
fprintf(dat, "cpd = 0.02; cphb= 0.75; cpho= 0.25; cpa = atan(256*cpd/2);\n");
fprintf(dat, "for idx = 0:255\n");
fprintf(dat, "   t = 0.5*(atan(cpd*(idx-256/2))/cpa+1);\n");
fprintf(dat, "   mag_blue(idx+1, 3) = (2-4*cphb)*t*t+(4*cphb-1)*t;\n");
fprintf(dat, "   mag_blue(idx+1, 1) = (2-4*cpho)*t*t+(4*cpho-1)*t;\n");
fprintf(dat, "   mag_blue(idx+1, 2) = mag_blue(idx+1, 1);\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% A few style variables ...\n");
fprintf(dat, "global col_fbg col_txt col_abg col_atx col_but col_btx fntsz;\n");
fprintf(dat, "col_fbg = [0.77 0.80 0.87]; col_txt = [0.00 0.00 0.00];\n");
fprintf(dat, "col_abg = [1.00 1.00 1.00]; col_atx = [0.00 0.00 0.00];\n");
fprintf(dat, "col_but = [0.00 0.47 0.48]; col_btx = [1.00 1.00 1.00];\n");
fprintf(dat, "fntsz = 10;\n");
fprintf(dat, "\n");
fprintf(dat, "%% The viewer window\n");
fprintf(dat, "global fgr;\n");
fprintf(dat, "fgr = figure('NumberTitle', 'off', 'Name', '%s', 'Color', col_fbg, ...\n", name);
fprintf(dat, "             'Units', 'pixels', 'Menubar', 'none');\n");
fprintf(dat, "\n");
fprintf(dat, "%% Adjust figure extension, respect aspect ratio of the data as far as possible,\n");
fprintf(dat, "%% Select figure position\n");
fprintf(dat, "dx = xend-xbeg; dy = yend-ybeg;\n");
fprintf(dat, "if dx<0.25*dy dx = 0.25*dy; end\n");
fprintf(dat, "if dy<0.25*dx dy = 0.25*dx; end\n");
fprintf(dat, "axw = round(sqrt(144400*dy/dx)); axh = round(axw*dx/dy);\n");
fprintf(dat, "global csf_xaw csf_yaw;\n");
fprintf(dat, "csf_yaw = axw; csf_xaw = axh;\n");
fprintf(dat, "fgw = axw+80+180; fgh = axh+2*40;\n");
fprintf(dat, "axw = axw/fgw; axh = axh/fgh;\n");
fprintf(dat, "axpx = 60/fgw; axpy = 40/fgh;\n");
fprintf(dat, "cmpx = axpx+axw+7/fgw; cmpw = 13/fgw;\n");
fprintf(dat, "bupx = cmpx+cmpw+45/fgw; buw = 1.0-10/fgw-bupx; busep = 2/fgw;\n");
fprintf(dat, "set(0, 'Units', 'pixels');\n");
fprintf(dat, "scrsz = get(0, 'ScreenSize');\n");
fprintf(dat, "if fgw > scrsz(3)*0.8 \n");
fprintf(dat, "   s = scrsz(3)*0.8/fgw; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "if fgh > scrsz(4)*0.8 \n");
fprintf(dat, "   s = scrsz(4)*0.8/fgh; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "fgpx = round(scrsz(1)+scrsz(3)*0.05); fgpy = round(scrsz(2)+scrsz(4)*0.95)-fgh;\n");
fprintf(dat, "set(fgr, 'Position', [fgpx fgpy fgw fgh]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Plot container: global state \n");
fprintf(dat, "global pl_fld pl_fldre pl_fldim pl_tp pl_min pl_max pl_amin pl_amax;\n");
fprintf(dat, "global pl_cp pl_img pl_afo pl_tih an_nump;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Initial plot: basic field component \n");
fprintf(dat, "if strcmp(pol, 'TE')\n");
fprintf(dat, "   pl_tp = 0.0; pl_cp = 'E_y'; pl_afo = 'f';\n");
fprintf(dat, "   pl_fldre = Ey_re; pl_fldim = Ey_im; pl_min = -max_E; pl_max =  max_E;\n");
fprintf(dat, "elseif strcmp(pol, 'TM')\n");
fprintf(dat, "   pl_tp = 0.0; pl_cp = 'H_y'; pl_afo = 'f';\n");
fprintf(dat, "   pl_fldre = Hy_re; pl_fldim = Hy_im; pl_min = -max_H; pl_max =  max_H;\n");
fprintf(dat, "else\n");
fprintf(dat, "   pl_tp = 0.0; pl_cp = 'E_z'; pl_afo = 'f';\n");
fprintf(dat, "   pl_fldre = Ez_re; pl_fldim = Ez_im; pl_min = -max_E; pl_max =  max_E;\n");
fprintf(dat, "end\n");
fprintf(dat, "an_nump = 10;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Refine the plot window\n");
fprintf(dat, "global axs;\n");
fprintf(dat, "axs = axes;\n");
fprintf(dat, "pl_fld = pl_fldre*cos(2*pi*pl_tp)-pl_fldim*sin(2*pi*pl_tp);\n");
fprintf(dat, "pl_amin = min(min(pl_fld)); pl_amax = max(max(pl_fld));\n");
fprintf(dat, "pl_img = imagesc(y, x, pl_fld, [pl_min, pl_max]);\n");
fprintf(dat, "set(axs, 'Position', [axpx axpy axw axh], 'Units', 'normalized', ...\n");
fprintf(dat, "         'YDir', 'normal', 'ButtonDownFcn', {@clbf_crossh}, ...\n");
fprintf(dat, "         'XLim', [ybeg, yend], 'YLim', [xbeg, xend], ...\n");
fprintf(dat, "         'TickDir', 'in', 'FontSize', fntsz);\n");
fprintf(dat, "xlabel(['y [' char(956) 'm]'], 'FontSize', fntsz); ylabel(['x [' char(956) 'm]'], 'FontSize', fntsz);\n");
fprintf(dat, "global tle;\n");
fprintf(dat, "tle = title(' ', 'FontSize', fntsz+2);\n");
fprintf(dat, "pl_tih = text(yend-(yend-ybeg)/30, xbeg+(xend-xbeg)/30, ...\n");
fprintf(dat, "              sprintf('t = %%#4.2f T', pl_tp), 'Color', 'w', ...\n");
fprintf(dat, "              'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "              'VerticalAlignment', 'bottom', 'FontSize', fntsz);\n");
fprintf(dat, "hold on; box on; \n");
fprintf(dat, "\n");
fprintf(dat, "%% The initial plot\n");
fprintf(dat, "set(pl_img, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "for idx=1:nhd\n");
fprintf(dat, "   plot([hdbeg(idx), hdend(idx)], [hdpos(idx), hdpos(idx)], ...\n");
fprintf(dat, "        'w-', 'LineWidth', 1.0, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "end\n");
fprintf(dat, "for idx=1:nvd\n");
fprintf(dat, "   plot([vdpos(idx), vdpos(idx)], [vdbeg(idx), vdend(idx)], ...\n");
fprintf(dat, "        'w-', 'LineWidth', 1.0, 'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "end\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Crosshair \n");
fprintf(dat, "global chh chv chx chy;\n");
fprintf(dat, "chx = xend; chy = yend;\n");
fprintf(dat, "axes(axs);\n");
fprintf(dat, "chh = plot([ybeg yend], [chx chx], ...\n");
fprintf(dat, "           'LineStyle', '-.', 'LineWidth', 0.5, 'Color', [0.5 0.5 0.5], ...\n");
fprintf(dat, "           'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "chv = plot([chy chy], [xbeg xend], ...\n");
fprintf(dat, "           'LineStyle', '-.', 'LineWidth', 0.5, 'Color', [0.5 0.5 0.5], ...\n");
fprintf(dat, "           'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "\n");
fprintf(dat, "%% The colormap axis\n");
fprintf(dat, "global clb clb_img;\n");
fprintf(dat, "clb = axes;\n");
fprintf(dat, "rge = linspace(pl_min, pl_max, 64);\n");
fprintf(dat, "clb_img = imagesc([0 1], rge, [rge',rge'], [pl_min, pl_max]);\n");
fprintf(dat, "set(clb, 'Position', [cmpx axpy cmpw axh], 'Units', 'normalized', ...\n");
fprintf(dat, "         'YDir', 'normal', 'ButtonDownFcn', {@clbf_contour}, ...\n");
fprintf(dat, "         'XLim', [0, 1], 'YLim', [pl_min, pl_max], 'XTick', [], ...\n");
fprintf(dat, "         'YAxisLocation', 'right', 'TickDir', 'in', 'FontSize', fntsz);\n");
fprintf(dat, "set(clb_img, 'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "\n");
fprintf(dat, "%% The contour positioner\n");
fprintf(dat, "global cont_ah cont_minh cont_maxh cont_hdl cont_num; \n");
fprintf(dat, "lv = pl_min;\n");
fprintf(dat, "cont_ah = plot([0 1], [lv lv], ...\n");
fprintf(dat, "               'LineStyle', '-', 'LineWidth', 0.5, 'Color', [0.5 0.5 0.5], ...\n");
fprintf(dat, "               'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "cont_minh = plot([0 1], [pl_amin pl_amin], ...\n");
fprintf(dat, "                 'LineStyle', '-', 'LineWidth', 0.5, 'Color', [1.0 1.0 1.0], ...\n");
fprintf(dat, "                 'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "cont_maxh = plot([0 1], [pl_amax pl_amax], ...\n");
fprintf(dat, "                 'LineStyle', '-', 'LineWidth', 0.5, 'Color', [1.0 1.0 1.0], ...\n");
fprintf(dat, "                 'ButtonDownFcn', {@clbf_contour});\n");
fprintf(dat, "cont_num = 0;\n");
fprintf(dat, "\n");
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% The gui control objects \n");
fprintf(dat, "\n");
fprintf(dat, "uc_dummy = uicontrol('Style', 'text', 'String', 'M', 'Units', 'normalized', ...\n");
fprintf(dat, "                     'FontSize', fntsz, 'Visible', 'off');\n");
fprintf(dat, "dext = get(uc_dummy, 'Extent'); delete(uc_dummy);\n");
fprintf(dat, "uch0 = dext(4);\n");
fprintf(dat, "uch = (axh-37*busep)/11;\n");
fprintf(dat, "if uch0>uch uch0 = uch; end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotation: y-coordinate \n");
fprintf(dat, "uch = uch0; ucpy = axpy; ucpx = bupx;\n");
fprintf(dat, "uc_yl = uicontrol('Style', 'text', 'String', 'y/um', 'Units', 'normalized', ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "p = get(uc_yl, 'Extent'); ucw = p(3);\n");
fprintf(dat, "set(uc_yl, 'Position', [ucpx ucpy ucw uch]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Input field: y-coordinate \n");
fprintf(dat, "global uc_yf;\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = buw-busep-ucw-busep-(buw-2*busep)/3;\n");
fprintf(dat, "uc_yf = uicontrol('Style', 'edit', 'Callback', {@cbf_yfld}, ...\n");
fprintf(dat, "                  'String', num2str(chy), 'Min', 1, 'Max', 1, ...\n");
fprintf(dat, "                  'TooltipString', 'View & set y-position', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_abg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: vertical cross section\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_vc = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_crosssec, 'v'}, ...\n");
fprintf(dat, "                  'String', '|', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Vertical field section', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotation: x-coordinate \n");
fprintf(dat, "ucpy = ucpy+uch+busep; ucpx = bupx;\n");
fprintf(dat, "uc_xl = uicontrol('Style', 'text', 'String', 'x/um', 'Units', 'normalized', ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "p = get(uc_xl, 'Extent'); ucw = p(3);\n");
fprintf(dat, "set(uc_xl, 'Position', [ucpx ucpy ucw uch]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Input field: x-coordinate \n");
fprintf(dat, "global uc_xf;\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = buw-busep-ucw-busep-(buw-2*busep)/3;\n");
fprintf(dat, "uc_xf = uicontrol('Style', 'edit', 'Callback', {@cbf_xfld}, ...\n");
fprintf(dat, "                  'String', num2str(chx), 'Min', 1, 'Max', 1, ...\n");
fprintf(dat, "                  'TooltipString', 'View & set x-position', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_abg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: horizontal cross section\n");
fprintf(dat, "ucpx = ucpx+busep+ucw; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_hc = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_crosssec, 'h'}, ...\n");
fprintf(dat, "                  'String', '---', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Horizontal field section', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: refractive index profile\n");
fprintf(dat, "ucw = (buw-2*busep)/3; ucpx = bupx+2*ucw+2*busep; ucpy = ucpy+uch+busep;\n");
fprintf(dat, "uc_n  = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_npl}, ...\n");
fprintf(dat, "                  'String', 'n', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Refractive index profile', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%%\n");
fprintf(dat, "ucpx = bupx; ucpy = ucpy+uch+10*busep; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "%% Button: Sz-profile \n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_sz = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_szpl}, ...\n");
fprintf(dat, "                  'String', 'Sz', 'FontSize', fntsz, ...\n");
fprintf(dat, "            'TooltipString', 'Time-averaged Poynting vector, z-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: w-profile \n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_w  = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_wpl}, ...\n");
fprintf(dat, "                  'String', 'w', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Energy density of the optical field', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: field-profile \n");
fprintf(dat, "ucpx = bupx; ucpy = ucpy+uch+busep; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_f  = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_fldpl}, ...\n");
fprintf(dat, "                  'String', 'F', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'The physical field, time snapshot', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: profile of the absolute field value\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_fm = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_fldmpl}, ...\n");
fprintf(dat, "                  'String', '|F|', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Absolute value of the complex field', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: profile of the absolute field square \n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_fq = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_fldqpl}, ...\n");
fprintf(dat, "                  'String', '|F|^2', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Absolute square of the complex field', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotations, slider time position\n");
fprintf(dat, "ucpy = ucpy+uch+busep; uch = uch0; ucw = buw/3; ucpx = bupx;\n");
fprintf(dat, "uc_tll = uicontrol('Style', 'text', 'String', 't=0 ...', ...\n");
fprintf(dat, "                   'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                   'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                   'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "ucpx = ucpx+ucw;\n");
fprintf(dat, "uc_tlc = uicontrol('Style', 'text', 'String', 'T/2', ...\n");
fprintf(dat, "                   'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                   'FontSize', fntsz, 'HorizontalAlignment', 'center', ...\n");
fprintf(dat, "                   'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "ucpx = ucpx+ucw;\n");
fprintf(dat, "uc_tlc = uicontrol('Style', 'text', 'String', '... T ', ...\n");
fprintf(dat, "                   'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                   'FontSize', fntsz, 'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "                   'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Slider: field snapshot time \n");
fprintf(dat, "global uc_tp;\n");
fprintf(dat, "ucpy = ucpy+uch+busep; uch = uch0*3/4; ucw = buw; ucpx = bupx;\n");
fprintf(dat, "uc_tp = uicontrol('Style', 'slider', 'Callback', {@cbf_timep}, ...\n");
fprintf(dat, "                  'Max', 1.0, 'Min', 0.0, 'Value', 0.0, ...\n");
fprintf(dat, "      'TooltipString', 'Physical field, snapshot time within one period T', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: play reverse animation\n");
fprintf(dat, "ucpy = ucpy+uch+busep; uch = uch0; ucpx = bupx; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_anbw = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_play, 'b'}, ...\n");
fprintf(dat, "                    'String', '<<', 'FontSize', fntsz, ...\n");
fprintf(dat, "                    'TooltipString', 'Backward wave propagation', ...\n");
fprintf(dat, "                    'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                    'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Input field: number of time periods shown in the animation\n");
fprintf(dat, "global uc_annp; ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_annp = uicontrol('Style', 'edit', 'Callback', {@cbf_annp}, ...\n");
fprintf(dat, "                    'String', num2str(an_nump), 'Min', 1, 'Max', 1, ...\n");
fprintf(dat, "                    'TooltipString', 'Runtime/T of the animation', ...\n");
fprintf(dat, "                    'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                    'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                    'BackgroundColor', col_abg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: play field animation\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_anfw = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_play, 'f'}, ...\n");
fprintf(dat, "                    'String', '>>', 'FontSize', fntsz, ...\n");
fprintf(dat, "                    'TooltipString', 'Forward wave propagation', ...\n");
fprintf(dat, "                    'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                    'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Hx-component\n");
fprintf(dat, "ucpy = ucpy+uch+10*busep; uch = uch0; ucpx = bupx; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_hx = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_hxpl}, ...\n");
fprintf(dat, "                  'String', 'Hx', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Magnetic field, x-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Hy-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_hy = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_hypl}, ...\n");
fprintf(dat, "                  'String', 'Hy', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Magnetic field, y-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Hz-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_hz = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_hzpl}, ...\n");
fprintf(dat, "                  'String', 'Hz', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Magnetic field, z-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Ex-component\n");
fprintf(dat, "ucpx = bupx; ucpy = ucpy+uch+busep; ucw = (buw-2*busep)/3;\n");
fprintf(dat, "uc_ex = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_expl}, ...\n");
fprintf(dat, "                  'String', 'Ex', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Electric field, x-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Ey-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_ey = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_eypl}, ...\n");
fprintf(dat, "                  'String', 'Ey', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Electric field, y-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Button: Ey-component\n");
fprintf(dat, "ucpx = ucpx+busep+ucw;\n");
fprintf(dat, "uc_ez = uicontrol('Style', 'pushbutton', 'Callback', {@cbf_ezpl}, ...\n");
fprintf(dat, "                  'String', 'Ez', 'FontSize', fntsz, ...\n");
fprintf(dat, "                  'TooltipString', 'Electric field, z-component', ...\n");
fprintf(dat, "                  'Units', 'normalized', 'Position', [ucpx ucpy ucw uch], ...\n");
fprintf(dat, "                  'BackgroundColor', col_but, 'ForegroundColor', col_btx);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Annotation: polarization and wavelength\n");
fprintf(dat, "uch = uch0; ucpy = axpy+axh-uch; ucpx = bupx;\n");
fprintf(dat, "uc_an = uicontrol('Style', 'text', 'Units', 'normalized', ...\n");
fprintf(dat, "                  'String', [pol ', ' num2str(lambda) ' um'], ...\n");
fprintf(dat, "                  'FontSize', fntsz, 'HorizontalAlignment', 'left', ...\n");
fprintf(dat, "                  'BackgroundColor', col_fbg, 'ForegroundColor', col_txt);\n");
fprintf(dat, "p = get(uc_an, 'Extent'); ucw = p(3);\n");
fprintf(dat, "set(uc_an, 'Position', [ucpx ucpy ucw uch]);\n");
fprintf(dat, "\n");
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% Callback functions for the gui objects \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: input y-coordinate \n");
fprintf(dat, "function cbf_yfld(obj, eventdata)\n");
fprintf(dat, "global uc_yf chy chv ybeg yend;\n");
fprintf(dat, "chy = str2double(get(uc_yf, 'String'));\n");
fprintf(dat, "if isnan(chy) chy = 0.5*(ybeg+yend); end\n");
fprintf(dat, "if chy<ybeg chy = ybeg; end\n");
fprintf(dat, "if chy>yend chy = yend; end\n");
fprintf(dat, "set(chv, 'XData', [chy chy]); \n");
fprintf(dat, "set(uc_yf, 'String', num2str(chy)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: input x-coordinate \n");
fprintf(dat, "function cbf_xfld(obj, eventdata)\n");
fprintf(dat, "global uc_xf chx chh xbeg xend;\n");
fprintf(dat, "chx = str2double(get(uc_xf, 'String'));\n");
fprintf(dat, "if isnan(chx) chx = 0.5*(xbeg+xend); end\n");
fprintf(dat, "if chx<xbeg chx = xbeg; end\n");
fprintf(dat, "if chx>xend chx = xend; end\n");
fprintf(dat, "set(chh, 'YData', [chx chx]); \n");
fprintf(dat, "set(uc_xf, 'String', num2str(chx)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Sz \n");
fprintf(dat, "function cbf_szpl(obj, eventdata)\n");
fprintf(dat, "global Sz max_S pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Sz; pl_fldim = Sz*0; pl_cp = 'S_z';\n");
fprintf(dat, "pl_min = -max_S; pl_max = max_S;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select w \n");
fprintf(dat, "function cbf_wpl(obj, eventdata)\n");
fprintf(dat, "global W max_W pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = W; pl_fldim = W*0; pl_cp = 'w'; pl_min = 0.0; pl_max = max_W;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select physical field \n");
fprintf(dat, "function cbf_fldpl(obj, eventdata)\n");
fprintf(dat, "global pl_afo;\n");
fprintf(dat, "pl_afo = 'f'; \n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select absolute field value\n");
fprintf(dat, "function cbf_fldmpl(obj, eventdata)\n");
fprintf(dat, "global pl_afo;\n");
fprintf(dat, "pl_afo = 'm'; \n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select absolute field square \n");
fprintf(dat, "function cbf_fldqpl(obj, eventdata)\n");
fprintf(dat, "global pl_afo;\n");
fprintf(dat, "pl_afo = 'q'; \n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select adjust snapshot time \n");
fprintf(dat, "function cbf_timep(obj, eventdata)\n");
fprintf(dat, "global uc_tp pl_tp;\n");
fprintf(dat, "pl_tp = get(uc_tp, 'Value');\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select animation runtime \n");
fprintf(dat, "function cbf_annp(obj, eventdata)\n");
fprintf(dat, "global uc_annp an_nump;\n");
fprintf(dat, "np = str2double(get(uc_annp, 'String'));\n");
fprintf(dat, "if isnan(np) np = 10; end\n");
fprintf(dat, "if np<1  np = 1;  end\n");
fprintf(dat, "if np>99 np = 99; end\n");
fprintf(dat, "an_nump = round(np);\n");
fprintf(dat, "set(uc_annp, 'String', num2str(an_nump)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Hx \n");
fprintf(dat, "function cbf_hxpl(obj, eventdata)\n");
fprintf(dat, "global Hx_re Hx_im max_H pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Hx_re; pl_fldim = Hx_im; pl_cp = 'H_x';\n");
fprintf(dat, "pl_min = -max_H; pl_max = max_H;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Hy \n");
fprintf(dat, "function cbf_hypl(obj, eventdata)\n");
fprintf(dat, "global Hy_re Hy_im max_H pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Hy_re; pl_fldim = Hy_im; pl_cp = 'H_y';\n");
fprintf(dat, "pl_min = -max_H; pl_max = max_H;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Hz \n");
fprintf(dat, "function cbf_hzpl(obj, eventdata)\n");
fprintf(dat, "global Hz_re Hz_im max_H pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Hz_re; pl_fldim = Hz_im; pl_cp = 'H_z';\n");
fprintf(dat, "pl_min = -max_H; pl_max = max_H;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Ex \n");
fprintf(dat, "function cbf_expl(obj, eventdata)\n");
fprintf(dat, "global Ex_re Ex_im max_E pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Ex_re; pl_fldim = Ex_im; pl_cp = 'E_x';\n");
fprintf(dat, "pl_min = -max_E; pl_max = max_E;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Ey \n");
fprintf(dat, "function cbf_eypl(obj, eventdata)\n");
fprintf(dat, "global Ey_re Ey_im max_E pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Ey_re; pl_fldim = Ey_im; pl_cp = 'E_y';\n");
fprintf(dat, "pl_min = -max_E; pl_max = max_E;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select Ez \n");
fprintf(dat, "function cbf_ezpl(obj, eventdata)\n");
fprintf(dat, "global Ez_re Ez_im max_E pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = Ez_re; pl_fldim = Ez_im; pl_cp = 'E_z';\n");
fprintf(dat, "pl_min = -max_E; pl_max = max_E;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: select n\n");
fprintf(dat, "function cbf_npl(obj, eventdata)\n");
fprintf(dat, "global n min_n max_n pl_fldre pl_fldim pl_min pl_max pl_cp;\n");
fprintf(dat, "pl_fldre = n; pl_fldim = n*0.0; pl_cp = 'n'; pl_min = min_n; pl_max = max_n;\n");
fprintf(dat, "if pl_min>=pl_max pl_min = pl_min-0.5; pl_max = pl_min+1.0; end;\n");
fprintf(dat, "d = (pl_max-pl_min)/3; pl_min = pl_min-d; pl_max = pl_max+d;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% -------------------------------------------------------------------------\n");
fprintf(dat, "%% The actual plotting routines, partly callbacks for gui objects\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: crosshair positioning \n");
fprintf(dat, "function clbf_crossh(obj, eventdata)\n");
fprintf(dat, "global fgr axs chh chv chx chy;\n");
fprintf(dat, "mpos = get(0, 'PointerLocation'); \n");
fprintf(dat, "mx=mpos(1); my=mpos(2);\n");
fprintf(dat, "fpos = get(fgr, 'Position'); fx=fpos(1); fy=fpos(2); fw=fpos(3); fh=fpos(4); \n");
fprintf(dat, "rx = (mx-fx)/fw; ry = (my-fy)/fh;\n");
fprintf(dat, "apos = get(axs, 'Position'); ax=apos(1); ay=apos(2); aw=apos(3); ah=apos(4); \n");
fprintf(dat, "extx = get(axs, 'Xlim'); x0=extx(1); x1=extx(2);\n");
fprintf(dat, "exty = get(axs, 'Ylim'); y0=exty(1); y1=exty(2);\n");
fprintf(dat, "px = x0+(rx-ax)/aw*(x1-x0); py = y0+(ry-ay)/ah*(y1-y0);\n");
fprintf(dat, "chx = py; chy = px;\n");
fprintf(dat, "global xbeg xend ybeg yend;\n");
fprintf(dat, "if chx<xbeg chx = xbeg; end\n");
fprintf(dat, "if chx>xend chx = xend; end\n");
fprintf(dat, "if chy<ybeg chy = ybeg; end\n");
fprintf(dat, "if chy>yend chy = yend; end\n");
fprintf(dat, "set(chh, 'YData', [chx chx]); set(chv, 'XData', [chy chy]); \n");
fprintf(dat, "global uc_xf uc_yf;\n");
fprintf(dat, "set(uc_xf, 'String', num2str(chx)); set(uc_yf, 'String', num2str(chy)); \n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: place contour line\n");
fprintf(dat, "function clbf_contour(obj, eventdata)\n");
fprintf(dat, "global fgr axs clb cont_ah cont_num cont_hdl;\n");
fprintf(dat, "for k=1:cont_num delete(cont_hdl(k)); end\n");
fprintf(dat, "mpos = get(0, 'PointerLocation'); my=mpos(2);\n");
fprintf(dat, "fpos = get(fgr, 'Position'); fy=fpos(2); fh=fpos(4); \n");
fprintf(dat, "ry = (my-fy)/fh;\n");
fprintf(dat, "apos = get(clb, 'Position'); ay=apos(2); ah=apos(4); \n");
fprintf(dat, "exty = get(clb, 'Ylim'); y0=exty(1); y1=exty(2);\n");
fprintf(dat, "lv = y0+(ry-ay)/ah*(y1-y0); \n");
fprintf(dat, "global pl_amin pl_amax;\n");
fprintf(dat, "if lv<pl_amin lv = pl_amin; end\n");
fprintf(dat, "if lv>pl_amax lv = pl_amax; end\n");
fprintf(dat, "set(cont_ah, 'YData', [lv lv]); \n");
fprintf(dat, "axes(axs); hold on;\n");
fprintf(dat, "global y x pl_fld;\n");
fprintf(dat, "cm = contourc(y, x, pl_fld, [lv, lv]); cl = size(cm, 2);\n");
fprintf(dat, "k = 1; cont_hdl = ones(1, cl); cont_num = 0;\n");
fprintf(dat, "while(k < cl)\n");
fprintf(dat, "      np = cm(2, k); py = cm(1, k+1:k+np); px = cm(2, k+1:k+np);\n");
fprintf(dat, "      cont_num = cont_num+1;\n");
fprintf(dat, "      cont_hdl(cont_num) = line(py, px, ...\n");
fprintf(dat, "                          'Color', 'b', 'LineStyle', '-', 'LineWidth', 0.5, ...\n");
fprintf(dat, "                          'ButtonDownFcn', {@clbf_crossh});\n");
fprintf(dat, "      k = k+np+1;\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: update the field plot \n");
fprintf(dat, "function updateplot\n");
fprintf(dat, "global fgr axs tle clb clb_img;\n");
fprintf(dat, "global pl_fld pl_fldre pl_fldim pl_tp pl_min pl_max pl_amin pl_amax;\n");
fprintf(dat, "global pl_cp pl_img pl_afo pl_tih;\n");
fprintf(dat, "global x y;\n");
fprintf(dat, "global cont_ah cont_minh cont_maxh cont_num cont_hdl;\n");
fprintf(dat, "global mag_blue lin_blue;\n");
fprintf(dat, "for k=1:cont_num delete(cont_hdl(k)); end\n");
fprintf(dat, "cont_num = 0;\n");
fprintf(dat, "switch pl_cp\n");
fprintf(dat, "case {'E_x', 'E_y', 'E_z', 'H_x', 'H_y', 'H_z'} \n");
fprintf(dat, "   if pl_afo == 'm' \n");
fprintf(dat, "      title = ['| ' pl_cp '(x, y) |']; \n");
fprintf(dat, "      pl_fld = sqrt(pl_fldre.*pl_fldre+pl_fldim.*pl_fldim);\n");
fprintf(dat, "      clm = lin_blue; minlv = 0.0; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   elseif pl_afo == 'q'\n");
fprintf(dat, "      title = ['| ' pl_cp '(x, y) |^2']; \n");
fprintf(dat, "      pl_fld = pl_fldre.*pl_fldre+pl_fldim.*pl_fldim;\n");
fprintf(dat, "      clm = lin_blue; minlv = 0.0; maxlv = pl_max*pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   else\n");
fprintf(dat, "      title = ['Re ' pl_cp '(x, y) e^{i2' char(960) 't/T}'];\n");
fprintf(dat, "      pl_fld = pl_fldre*cos(2*pi*pl_tp)-pl_fldim*sin(2*pi*pl_tp);\n");
fprintf(dat, "      clm = mag_blue; minlv = -pl_max; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'on', 'String', sprintf('t = %%#4.2f T', pl_tp));\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'S_x', 'S_y', 'S_z'}\n");
fprintf(dat, "   if pl_afo == 'm' || pl_afo == 'q'\n");
fprintf(dat, "      title = ['| ' pl_cp '(x, y) |']; \n");
fprintf(dat, "      pl_fld = sqrt(pl_fldre.*pl_fldre);\n");
fprintf(dat, "      clm = lin_blue; minlv = 0.0; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   else\n");
fprintf(dat, "      title = [pl_cp '(x, y)'];\n");
fprintf(dat, "      pl_fld = pl_fldre;\n");
fprintf(dat, "      clm = mag_blue; minlv = -pl_max; maxlv = pl_max;\n");
fprintf(dat, "      set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'w'}\n");
fprintf(dat, "   title = 'w(x, y)';\n");
fprintf(dat, "   pl_fld = pl_fldre;\n");
fprintf(dat, "   clm = lin_blue; minlv = 0.0; maxlv = pl_max;\n");
fprintf(dat, "   set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "case {'n'}\n");
fprintf(dat, "   title = 'n(x, y)';\n");
fprintf(dat, "   pl_fld = pl_fldre;\n");
fprintf(dat, "   clm = lin_blue; minlv = pl_min; maxlv = pl_max;\n");
fprintf(dat, "   set(pl_tih, 'Visible', 'off');\n");
fprintf(dat, "end\n");
fprintf(dat, "if minlv == maxlv\n");
fprintf(dat, "   minlv = minlv-1.0\n");
fprintf(dat, "   maxlv = minlv+2.0\n");
fprintf(dat, "end\n");
fprintf(dat, "pl_amin = min(min(pl_fld)); pl_amax = max(max(pl_fld));\n");
fprintf(dat, "set(fgr, 'Colormap', clm);\n");
fprintf(dat, "set(clb_img, 'CData', linspace(minlv, maxlv, 64)', ...\n");
fprintf(dat, "             'YData', linspace(minlv, maxlv, 64)');\n");
fprintf(dat, "set(clb, 'YLim', [minlv, maxlv], 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "set(cont_ah, 'YData', [minlv minlv]); \n");
fprintf(dat, "set(cont_minh, 'YData', [pl_amin pl_amin]); \n");
fprintf(dat, "set(cont_maxh, 'YData', [pl_amax pl_amax]); \n");
fprintf(dat, "set(tle, 'String', title);\n");
fprintf(dat, "set(axs, 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "set(pl_img, 'CData', pl_fld);\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: field animation \n");
fprintf(dat, "function cbf_play(obj, eventdata, dr) \n");
fprintf(dat, "global fgr axs tle clb clb_img uc_tp mag_blue lin_blue;\n");
fprintf(dat, "global pl_fld pl_fldre pl_fldim pl_tp pl_min pl_max x y;\n");
fprintf(dat, "global pl_cp pl_img pl_afo pl_tih an_nump;\n");
fprintf(dat, "global cont_ah cont_num cont_hdl;\n");
fprintf(dat, "for k=1:cont_num delete(cont_hdl(k)); end\n");
fprintf(dat, "cont_num = 0;\n");
fprintf(dat, "switch pl_cp\n");
fprintf(dat, "case {'E_x', 'E_y', 'E_z', 'H_x', 'H_y', 'H_z'} \n");
fprintf(dat, "   if pl_afo == 'm' | pl_afo =='q' return; end\n");
fprintf(dat, "case {'S_x', 'S_y', 'S_z', 'w', 'n'}\n");
fprintf(dat, "   return;\n");
fprintf(dat, "end\n");
fprintf(dat, "clm = mag_blue;\n");
fprintf(dat, "minlv = -pl_max; maxlv = pl_max;\n");
fprintf(dat, "set(fgr, 'Colormap', clm);\n");
fprintf(dat, "set(clb_img, 'CData', linspace(minlv, maxlv, 64)', ...\n");
fprintf(dat, "             'YData', linspace(minlv, maxlv, 64)');\n");
fprintf(dat, "set(clb, 'YLim', [minlv, maxlv], 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "set(cont_ah, 'YData', [minlv minlv]); \n");
fprintf(dat, "set(axs, 'CLim', [minlv, maxlv]);\n");
fprintf(dat, "title = ['Re ' pl_cp '(x, y) e^{i2' char(960) 't/T}'];\n");
fprintf(dat, "set(tle, 'String', title);\n");
fprintf(dat, "set(pl_tih, 'Visible', 'on', 'String', sprintf('t = %%#4.2f T', pl_tp));\n");
fprintf(dat, "axes(axs);\n");
fprintf(dat, "msp = moviein(30);\n");
fprintf(dat, "for ti=0:29\n");
fprintf(dat, "   if dr == 'b'\n");
fprintf(dat, "      t = pl_tp-ti/30; if t<0 t = 1+t; end\n");
fprintf(dat, "   else \n");
fprintf(dat, "      t = pl_tp+ti/30; if t>1 t = t-1; end\n");
fprintf(dat, "   end\n");
fprintf(dat, "   pl_fld = pl_fldre*cos(2*pi*t)-pl_fldim*sin(2*pi*t);\n");
fprintf(dat, "   set(pl_img, 'CData', pl_fld);\n");
fprintf(dat, "   set(pl_tih, 'String', sprintf('t = %%#4.2f T', t));\n");
fprintf(dat, "   set(uc_tp, 'Value', t);\n");
fprintf(dat, "   msp(:,ti+1) = getframe;\n");
fprintf(dat, "end\n");
fprintf(dat, "set(uc_tp, 'Value', pl_tp);\n");
fprintf(dat, "movie(msp, an_nump-1);\n");
fprintf(dat, "title = ['Re ' pl_cp '(x, y) e^{i2' char(960) 't/T}'];\n");
fprintf(dat, "set(tle, 'String', title);\n");
fprintf(dat, "pl_fld = pl_fldre*cos(2*pi*(pl_tp+t))-pl_fldim*sin(2*pi*(pl_tp+t));\n");
fprintf(dat, "set(pl_img, 'CData', pl_fld);\n");
fprintf(dat, "set(pl_tih, 'String', sprintf('t = %%#4.2f T', pl_tp));\n");
fprintf(dat, "clear msp;\n");
fprintf(dat, "updateplot;\n");
fprintf(dat, "\n");
fprintf(dat, "%% Callback: cross section figure\n");
fprintf(dat, "function cbf_crosssec(obj, eventdata, dr)\n");
fprintf(dat, "global csf_xaw csf_yaw;\n");
fprintf(dat, "global col_fbg col_txt col_abg col_atx fntsz fgr;\n");
fprintf(dat, "global x y xbeg xend ybeg yend numx numy;\n");
fprintf(dat, "global chx chy pl_fld pl_cp pl_afo pl_tp pl_min pl_max; \n");
fprintf(dat, "csf = figure('NumberTitle', 'off', 'Name', '%s, field section', ...\n", name);
fprintf(dat, "             'Color', col_fbg, 'Units', 'pixels', 'Menubar', 'none');\n");
fprintf(dat, "if dr == 'h' axw = csf_yaw; else axw = csf_xaw; end\n");
fprintf(dat, "axh = 200; fgw = axw+60+20; fgh = axh+40+20;\n");
fprintf(dat, "axw = axw/fgw; axh = axh/fgh; axpx = 60/fgw; axpy = 40/fgh;\n");
fprintf(dat, "scrsz = get(0, 'ScreenSize');\n");
fprintf(dat, "if fgw > scrsz(3)*0.6 \n");
fprintf(dat, "   s = scrsz(3)*0.6/fgw; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "if fgh > scrsz(4)*0.6 \n");
fprintf(dat, "   s = scrsz(4)*0.6/fgh; fgw = fgw*s; fgh = fgh*s;\n");
fprintf(dat, "end\n");
fprintf(dat, "fgpx = round(scrsz(1)+scrsz(3)*0.1); fgpy = round(scrsz(2)+scrsz(4)*0.90)-fgh;\n");
fprintf(dat, "p = get(fgr, 'Position');\n");
fprintf(dat, "if p(1)+p(3)+fgw+10 <= scrsz(1)+scrsz(3) \n");
fprintf(dat, "   fgpx = p(1)+p(3)+10; fgpy = p(2);\n");
fprintf(dat, "end\n");
fprintf(dat, "set(csf, 'Position', [fgpx fgpy fgw fgh]);\n");
fprintf(dat, "csf_axs = axes;\n");
fprintf(dat, "axes(csf_axs); hold on; box on;\n");
fprintf(dat, "set(csf_axs, 'Position', [axpx axpy axw axh], 'Units', 'normalized', ...\n");
fprintf(dat, "             'YDir', 'normal', 'TickDir', 'out', 'FontSize', fntsz);\n");
fprintf(dat, "if dr == 'h'\n");
fprintf(dat, "   xlabel(['y [' char(956) 'm]'], 'FontSize', fntsz);\n");
fprintf(dat, "   cbeg = ybeg; cend = yend; cod = y;\n");
// fprintf(dat, "   poss = sprintf('x = %%.5g \\\\mum', chx);\n");
fprintf(dat, "   poss = sprintf(['x = %%.5g ' char(956) 'm'], chx);\n");
fprintf(dat, "else \n");
fprintf(dat, "   xlabel(['x [' char(956) 'm]'], 'FontSize', fntsz);\n");
fprintf(dat, "   cbeg = xbeg; cend = xend; cod = x;\n");
// fprintf(dat, "   poss = sprintf('y = %%.5g \\\\mum', chy);\n");
fprintf(dat, "   poss = sprintf(['y = %%.5g ' char(956) 'm'], chy);\n");
fprintf(dat, "end\n");
fprintf(dat, "set(csf_axs, 'XLim', [cbeg, cend]);\n");
fprintf(dat, "coc = linspace(cbeg, cend, 1000);\n");
fprintf(dat, "if dr == 'h'\n");
fprintf(dat, "   fldicoc = interp2(y, x, pl_fld, coc, linspace(chx, chx, 1000), 'linear');\n");
fprintf(dat, "   fldicod = interp2(y, x, pl_fld, cod, linspace(chx, chx, length(y)));\n");
fprintf(dat, "else\n");
fprintf(dat, "   fldicoc = interp2(y, x, pl_fld, linspace(chy, chy, 1000), coc, 'linear');\n");
fprintf(dat, "   fldicod = interp2(y, x, pl_fld, linspace(chy, chy, length(x)), cod);\n");
fprintf(dat, "end\n");
fprintf(dat, "switch pl_cp\n");
fprintf(dat, "case {'E_x', 'E_y', 'E_z', 'H_x', 'H_y', 'H_z'} \n");
fprintf(dat, "   if pl_afo == 'm' \n");
fprintf(dat, "      ylabel(['| ' pl_cp ' |'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = 0.0; fend = pl_max;\n");
fprintf(dat, "   elseif pl_afo == 'q'\n");
fprintf(dat, "      ylabel(['| ' pl_cp ' |^2'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = 0.0; fend = pl_max*pl_max;\n");
fprintf(dat, "   else\n");
fprintf(dat, "      ylabel(['Re ' pl_cp ' e^{i2' char(960) 't/T}'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = -pl_max; fend = pl_max;\n");
fprintf(dat, "      poss = [poss, sprintf(',  t = %%#4.2f T', pl_tp)];\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'S_x', 'S_y', 'S_z'}\n");
fprintf(dat, "   if pl_afo == 'm' || pl_afo == 'q'\n");
fprintf(dat, "      ylabel(['| ' pl_cp ' |'], 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = 0.0; fend = pl_max;\n");
fprintf(dat, "   else\n");
fprintf(dat, "      ylabel(pl_cp, 'FontSize', fntsz);\n");
fprintf(dat, "      fbeg = -pl_max; fend = pl_max;\n");
fprintf(dat, "   end\n");
fprintf(dat, "case {'w'}\n");
fprintf(dat, "   ylabel('w', 'FontSize', fntsz);\n");
fprintf(dat, "   fbeg = 0.0; fend = pl_max;\n");
fprintf(dat, "case {'n'}\n");
fprintf(dat, "   ylabel('n', 'FontSize', fntsz);\n");
fprintf(dat, "   fbeg = pl_min; fend = pl_max;\n");
fprintf(dat, "end\n");
fprintf(dat, "global n nvd vdbeg vdend vdpos nhd hdbeg hdend hdpos min_n max_n;\n");
fprintf(dat, "if dr == 'h'\n");
fprintf(dat, "   cb = ones(nvd); ncb = 0;\n");
fprintf(dat, "   for d=1:nvd\n");
fprintf(dat, "      if vdbeg(d) <= chx && chx < vdend(d)\n");
fprintf(dat, "         ncb = ncb+1; cb(ncb) = vdpos(d);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for d=1:ncb+1\n");
fprintf(dat, "      if d==1     b0=cbeg; else b0=cb(d-1); end\n");
fprintf(dat, "      if d==ncb+1 b1=cend; else b1=cb(d);   end\n");
fprintf(dat, "      rip = 0.5*(b1+b0);\n");
fprintf(dat, "      if ybeg <= rip && rip <= yend\n");
fprintf(dat, "         ri = interp2(y, x, n, 0.5*(b1+b0), chx);\n");
fprintf(dat, "         pc = 1.0-(ri-min_n)/(max_n-min_n)*0.2;\n");
fprintf(dat, "         patch([b0 b0 b1 b1], [fend fbeg fbeg fend], [pc pc pc], ...\n");
fprintf(dat, "              'LineStyle', 'none');\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "else\n");
fprintf(dat, "   cb = ones(nhd); ncb = 0;\n");
fprintf(dat, "   for d=1:nhd\n");
fprintf(dat, "      if hdbeg(d) <= chy && chy < hdend(d)\n");
fprintf(dat, "         ncb = ncb+1; cb(ncb) = hdpos(d);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "   for d=1:ncb+1\n");
fprintf(dat, "      if d==1     b0=cbeg; else b0=cb(d-1); end\n");
fprintf(dat, "      if d==ncb+1 b1=cend; else b1=cb(d);   end\n");
fprintf(dat, "      rip = 0.5*(b1+b0);\n");
fprintf(dat, "      if xbeg <= rip && rip <= xend\n");
fprintf(dat, "         ri = interp2(y, x, n, chy, 0.5*(b1+b0));\n");
fprintf(dat, "         pc = 1.0-(ri-min_n)/(max_n-min_n)*0.2;\n");
fprintf(dat, "         patch([b0 b0 b1 b1], [fend fbeg fbeg fend], [pc pc pc], ...\n");
fprintf(dat, "              'LineStyle', 'none');\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "for d=1:ncb\n");
fprintf(dat, "   plot([cb(d), cb(d)], [fbeg fend], 'k-', 'LineWidth', 1);\n");
fprintf(dat, "end\n");
fprintf(dat, "text(cend-(cend-cbeg)/30, fend-(fend-fbeg)/30, ...\n");
fprintf(dat, "     poss, 'Color', 'k', 'HorizontalAlignment', 'right', ...\n");
fprintf(dat, "     'VerticalAlignment', 'top', 'FontSize', fntsz);\n");
fprintf(dat, "set(csf_axs, 'YLim', [fbeg, fend]);\n");
fprintf(dat, "plot(coc, fldicoc, 'b-', 'LineWidth', 2); \n");
fprintf(dat, "ms = 2; if length(cod) <= 60 ms = 4; end\n");
fprintf(dat, "plot(cod, fldicod, 'ko', 'MarkerSize', ms, 'LineWidth', 2); \n");
fprintf(dat, "plot([cbeg cbeg cend cend cbeg], [fend fbeg fbeg fend fend], 'k-');\n");
fprintf(dat, "\n");
}