/* 
 * WMM-Application:
 * square waveguide surrounded by air, fundamental hybrid mode
 */

/*
 * WMM
 * Wave-matching method for mode analysis of dielectric waveguides
 * Manfred Lohmeyer 
 * University of Osnabrueck, Department of Physics
 * Barbarastrasse 7, D-49069 Osnabrueck, Germany
 * (1999)
 */

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

/* waveguide parameters */
#define Wgpni 1.5   // core refractive index
#define Wgpno 1.0   // cover: air
#define Wgpa  1.0   // width of the rectangular core 
#define Wgpl  1.5   // vacuum wavelength 

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

	g.hx(0) = -Wgpa/2.0;
	g.hx(1) =  Wgpa/2.0;
	g.hy(0) = -Wgpa/2.0;
	g.hy(1) =  Wgpa/2.0;

	g.n(0,0) = Wgpno;
	g.n(0,1) = Wgpno;
	g.n(0,2) = Wgpno;
	g.n(1,0) = Wgpno;
	g.n(1,1) = Wgpni;
	g.n(1,2) = Wgpno;
	g.n(2,0) = Wgpno;
	g.n(2,1) = Wgpno;
	g.n(2,2) = Wgpno;

	g.lambda = Wgpl;

	return g;
}

/* 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   = 8;
	p.ini_alpha_max = 2.0;

	p.ini_steps = 10;
	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;
}

int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	WMM_Mode m;
	double d = 0.0001;
		
	Rect disp(-1.0, -1.0, 1.0, 1.0);

	Waveguide wg = wgdef();
	
	// find one of the fundamental vectorial modes
	WMM_findfundmode(wg, VEC, SYM, 0.0, 0.0, par, '-', '-', ma);
	m = ma(0);
	// save mode
	m.write_def('0', '0');

	// profile data for all six electromagnetic field components
	m.acmfile('0', '0', disp, 70, 70);

	// fields evaluated on various lines on the cross section
	int f;
	Fcomp fc = EY;
	char pltype;
	for(f=0; f<=5; ++f)
	{
		switch(f)
		{
		case 0: fc = EX; break;
		case 1: fc = EY; break;
		case 2: fc = EZ; break;
		case 3: fc = HX; break;
		case 4: fc = HY; break;
		case 5: fc = HZ; break;
		}
		if(f==0 || f== 3) pltype = 'L'; 
		else              pltype = 'V'; 
		// y - axis
	    m.secmfile(fc,'h','0',500,pltype,0.0,-1.0, 0.0, 1.0);
		// x - axis
	    m.secmfile(fc,'v','0',500,pltype,-1.0,0.0, 1.0, 0.0);
		// top core boundary, inside the core
	    m.secmfile(fc,'i','0',500,pltype,Wgpa/2.0-d,-1.0,Wgpa/2.0-d,1.0);
		// top core boundary, outside the core
	    m.secmfile(fc,'o','0',500,'V',   Wgpa/2.0+d,-1.0,Wgpa/2.0+d,1.0);
		// right core boundary, inside the core
	    m.secmfile(fc,'l','0',500,pltype,-1.0,Wgpa/2.0-d,1.0,Wgpa/2.0-d);
		// right core boundary, outside the core
	    m.secmfile(fc,'r','0',500,'V',   -1.0,Wgpa/2.0+d,1.0,Wgpa/2.0+d);
	}

	ma.clear();

	return 0;
}