/* 
 * WMM-Application:
 * a simple rib waveguide
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"wmminc.h"

/* waveguide parameters */
#define Wgpns 1.9  // substrate refractive index
#define Wgpnf 2.2  // film refractive index
#define Wgpnc 1.0  // cover: air
#define Wgpw  1.0  // rib width
#define Wgph  0.5  // rib height
#define Wgpt  0.5  // remaining film thickness
#define Wgpl  1.3  // vacuum wavelength 

/* waveguide definition */
Waveguide wgdef()
{
	Waveguide g(2, 1);

	g.hx(0) = 0.0;
	g.hx(1) = Wgpt;
	g.hx(2) = Wgph+Wgpt;
	g.hy(0) = -Wgpw/2.0;
	g.hy(1) =  Wgpw/2.0;

	g.n(0,0) = Wgpns;
	g.n(0,1) = Wgpns;
	g.n(0,2) = Wgpns;
	g.n(1,0) = Wgpnf;
	g.n(1,1) = Wgpnf;
	g.n(1,2) = Wgpnf;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = Wgpnf;
	g.n(2,2) = Wgpnc;
	g.n(3,0) = Wgpnc;
	g.n(3,1) = Wgpnc;
	g.n(3,2) = Wgpnc;

	g.lambda = Wgpl;

	return g;
}

/* WMM analysis parameters */
WMM_Parameters pardef()
{
	WMM_Parameters p;

	p.vform = HXHY;
	p.vnorm = NRMMH;
	p.ccomp = CCALL;

	p.ini_d_alpha   = 0.05;
	p.ini_N_alpha   = 15;
	p.ini_alpha_max = 2.0;

	p.ini_steps = 30;
	p.ref_num   = 5;
	p.ref_exp   = 4.0;
	p.ref_sdf   = 0.5;

	p.fin_d_alpha   = 0.01;
	p.fin_N_alpha   = 30;
	p.fin_alpha_max = 3.0;

	p.btol   = 1.0e-7;
	p.mshift = 1.0e-10;

	return p;
}

/* calculate the fundamental vectorial mode,
   save it to a file, 
   write profile data */
int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	WMM_Mode m;
	double h, v;
		
	// define the waveguide
	Waveguide wg = wgdef();

	// find its fundamental mode, save it
	WMM_findfundmode(wg, VEC, SYM, 0.0, 0.0, par, '-', '-', ma);
	m = ma(0);
	m.write_def('-', '-');

	// the second run should start here	
//	m.read_def(VEC, SYM, '-', '-');

	// display window
	Rect display(-Wgph, -1.5*Wgpw, (Wgph+Wgpt)*1.2, 1.5*Wgpw);

	// the dominant magnetic mode field component, various views  
	m.mfile(HX, ORG, display, 75, 115, '0', '0', 'S');
	m.mfile(HX, MOD, display, 75, 115, '0', '0', 'C');
	m.mfile(HX, SQR, display, 75, 115, '0', '0', 'I');

	// six components in one plot
	m.acmfile('0', '0', display, 70, 115);

	// some sections through the dominant electric field profile,  
	// note the discontinuity properties  
	h = Wgpt+Wgph/2.0;
	m.secmfile(EY, 'h', 'a', 500, 'L', h, display.y0, h, display.y1);
//	m.writesec(EY, 'h', 'a', 500, h, display.y0, h, display.y1);
	h = Wgpt/2.0;
	m.secmfile(EY, 'h', 'b', 500, 'V', h, display.y0, h, display.y1);
//	m.writesec(EY, 'h', 'b', 500, h, display.y0, h, display.y1);
	h = -Wgpt/2.0;
	m.secmfile(EY, 'h', 'c', 500, 'V', h, display.y0, h, display.y1);
//	m.writesec(EY, 'h', 'c', 500, h, display.y0, h, display.y1);
	h = Wgph+Wgpt+0.1;
	m.secmfile(EY, 'h', 'd', 500, 'V', h, display.y0, h, display.y1);
//	m.writesec(EY, 'h', 'd', 500, h, display.y0, h, display.y1);
	v = 0.0;
	m.secmfile(EY, 'v', 'a', 500, 'L', display.x0, v, display.x1, v);
//	m.writesec(EY, 'v', 'a', 500, display.x0, v, display.x1, v);
	v = -Wgpw;
	m.secmfile(EY, 'v', 'b', 500, 'V', display.x0, v, display.x1, v);
//	m.writesec(EY, 'v', 'b', 500, display.x0, v, display.x1, v);
	v = -Wgpw/2.0-(1.0e-8);
	m.secmfile(EY, 'v', 'c', 500, 'V', display.x0, v, display.x1, v);
//	m.writesec(EY, 'v', 'c', 500, display.x0, v, display.x1, v);
	v = -Wgpw/2.0+(1.0e-8);
	m.secmfile(EY, 'v', 'd', 500, 'V', display.x0, v, display.x1, v);
//	m.writesec(EY, 'v', 'd', 500, display.x0, v, display.x1, v);

	// the WMM logo :-)
	m.fancymfile(HX, display, 150, 200, '0', '0');

	return 0;
}