/* 
 * WMM-Application:
 * a series of rib waveguides with varying etching depth, a common benchmark
 * problem 
 */

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

/* common waveguide parameters */
#define Wgpw  3.0   // rib width 
#define WgpT  1.0   // film thickness before etching
#define Wgpns 3.40  // substrate refractive index 
#define Wgpnf 3.44  // film refractive index 
#define Wgpnc 1.0   // cover: air
#define Wgpl  1.15  // vacuum wavelength 

/* waveguide definition 
   t:  remaining film thickness outside the rib */
Waveguide wgdef(double t)
{
	Waveguide g(2, 1);

	g.hx(0) = 0.0;
	g.hx(1) = t;
	g.hx(2) = 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   = 20;
	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.005;
	p.fin_N_alpha   = 30;
	p.fin_alpha_max = 3.0;

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

	return p;
}

/* calculate fundamental semivectorial modes of both polarizations,
   save modes, write profile data */
int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	WMM_Mode m;
	Polarization pol;
	int p, i;	
	Fcomp fc;
	Waveguide wg;

	// display region 
	Rect display(-1.0, -3.0,  1.3, 3.0);
	
	// compute TE and TM modes
	for(p=0; p<=1; ++p)
	{
		if(p==0) pol = QTE; 
		else     pol = QTM; 
		
		// five different etching depths
		for(i=1; i<=9; i+=2)
		{
			// define the waveguide
			wg = wgdef(i*0.1);
			
			// find the fundamental mode, save it
			WMM_findfundmode(wg, pol, SYM, 0.0, 0.0, 
			                 par, '0'+i, '-', ma);
			m = ma(0); ma.clear();
			m.write_def('0'+i, '-');

			// read it again	
			// m.read_def(pol, SYM, '0'+i, '-');

			// make a contour plot of the mode profile  
			fc = principalcomp(pol);
			m.mfile(fc, SQR, display, 100, 150, '0'+i, '0', 'C');
		}
	}
	return 0;
}