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

/*
 * wmmvis.cpp
 * output to MATLAB .m -files, visualization  
 * helper routines for mode.cpp and wmmarr.cpp
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"waveg.h"
#include"maxweq.h"
#include"wmmvis.h"

/* some global settings */

/* output uses colour */
int Mlop_Colour=YES;
/* output files include print commands */
int Mlop_Print=NO;

/* waveguide geometry and field extrema */
void mlout_geo(FILE *dat, Waveguide 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, Waveguide 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, Waveguide iwg, Waveguide pwg, Waveguide 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 = iwg.xyplane.y1;
		if(pp <= pwg.ny) py = pwg.hy(pp);
		else             py = pwg.xyplane.y1;
		if(po <= owg.ny) oy = owg.hy(po);
		else             oy = owg.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;
}

/* matrix of field values */
void mlout_fld(FILE *dat, int npx, int npy, Fcomp cp, Dmatrix fld)
{
	int xi, yi;
	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(yi=0; yi<=npy-2; ++yi)
//			fprintf(dat, "%g ...\n", fld(xi, yi));
			fprintf(dat, "%g ", fld(xi, yi));
		fprintf(dat, "%g ;\n", fld(xi, npy-1));
	}
	fprintf(dat, "];\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(const char *name, FILE *dat, Fcomp cp, Afo foa)
{
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(['y [' char(956) 'm]'], 'FontSize', bfs*1.1);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', bfs*1.1);\n");
	switch(foa)
	{
		case MOD:  
fprintf(dat, "title('|%c_%c|', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case ORG: 
fprintf(dat, "title('%c_%c', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case SQR:  
fprintf(dat, "title('|%c_%c|^2', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case REP:  
fprintf(dat, "title('Re %c_%c', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case IMP:  
fprintf(dat, "title('Im %c_%c', 'FontSize', bfs*1.3);\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;
}

/* surface plot */
void mlout_surface(const char *name, FILE *dat, Fcomp cp, Afo foa)
{
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(['y [' char(956) 'm]'], 'FontSize', 17);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', 17);\n");
	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;
}

/* intensity image */
void mlout_image(const char *name, FILE *dat, Fcomp cp, Afo foa)
{
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, "%% uncomment to to 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));
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, "axis equal;\n");
fprintf(dat, "box on;\n");
fprintf(dat, "bfs = abs((xend-xbeg)/(yend-ybeg));\n");
fprintf(dat, "if bfs>1 bfs=1/bfs; end;\n");
fprintf(dat, "bfs=10+5*bfs;\n");
fprintf(dat, "axis([ybeg, yend, xbeg, xend]);\n");
fprintf(dat, "set(gca, 'FontSize', bfs, ...\n");
fprintf(dat, "	 'LineWidth', 1.5, ...\n");
fprintf(dat, "	 'Clipping', 'off', ...\n");
fprintf(dat, "	 'YDir', 'normal', ...\n");
fprintf(dat, "	 'TickDir', 'out');\n");
fprintf(dat, "xlabel(['y [' char(956) 'm]'], 'FontSize', bfs*1.1);\n");
fprintf(dat, "ylabel(['x [' char(956) 'm]'], 'FontSize', bfs*1.1);\n");
	switch(foa)
	{
		case MOD:  
fprintf(dat, "title('|%c_%c|', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case ORG: 
fprintf(dat, "title('%c_%c', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case SQR:  
fprintf(dat, "title('|%c_%c|^2', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case REP:  
fprintf(dat, "title('Re %c_%c', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
		case IMP:  
fprintf(dat, "title('Im %c_%c', 'FontSize', bfs*1.3);\n",
             fldchr(cp), cpchr(cp));
			break;
	}
fprintf(dat, "\n");
fprintf(dat, "%% Add sketch of permittivity boundaries\n");
fprintf(dat, "hold on;\n");
fprintf(dat, "lw = 2.5;\n");
	if(Mlop_Colour==YES)
fprintf(dat, "ls = 'k-';\n");
	else
fprintf(dat, "ls = 'w-';\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], ls, '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], ls, 'LineWidth', lw);\n");
fprintf(dat, "      end\n");
fprintf(dat, "   end\n");
fprintf(dat, "end\n");
fprintf(dat, "\n");
	return;
}

/* six component surface plot */  
void mlout_acmfile(const 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(const 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;
}

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;
}

/* fancy plot */
void mlout_fancy(const 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.8;\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");
	if(Mlop_Colour==YES)
	{
fprintf(dat, "%% Choose a 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, "%% Choose a colormap\n");
fprintf(dat, "gr = zeros(64, 3);\n");
fprintf(dat, "gr(:,1) = 0.7';\n");
fprintf(dat, "gr(:,2) = gr(:,1);\n");
fprintf(dat, "gr(:,3) = gr(:,1);\n");
fprintf(dat, "colormap(gr);\n");
	}
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, "\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, const char *name, const 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(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");
	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(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");
	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, const 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;	
}