/* 
 * WMM-Application:
 * three rib waveguides, a frequently used benchmark problem 
 */

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

/* common waveguide parameters */
#define Wgpnc 1.0  // cover: air
#define Wgpl  1.55 // vacuum wavelength 

/* waveguide definition 
   ns: substrate refractive index
   nf: film refractive index
   w:  rib width
   h:  rib height
   t:  remaining film thickness */
Waveguide wgdef(double ns, double nf, double w, double h, double t)
{
	Waveguide g(2, 1);

	g.hx(0) = 0.0;
	g.hx(1) = t;
	g.hx(2) = h+t;
	g.hy(0) = -w/2.0;
	g.hy(1) =  w/2.0;

	g.n(0,0) = ns;
	g.n(0,1) = ns;
	g.n(0,2) = ns;
	g.n(1,0) = nf;
	g.n(1,1) = nf;
	g.n(1,2) = nf;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = nf;
	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 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;

	// three waveguides and display regions 
	Waveguide wg[3];
	wg[0] = wgdef(3.340, 3.44, 2.0, 1.10, 0.20);
	wg[1] = wgdef(3.360, 3.44, 3.0, 0.10, 0.90);
	wg[2] = wgdef(3.435, 3.44, 4.0, 2.50, 3.50);
	Rect display[3];
	display[0] = Rect(-0.5, -1.7,  1.5, 1.7);
	display[1] = Rect(-1.0, -4.0,  1.5, 4.0);
	display[2] = Rect(-3.5, -10.0, 7.0, 10.0);
	
	// compute TE and TM modes
	for(p=0; p<=1; ++p)
	{
		if(p==0) pol = QTE; 
		else     pol = QTM; 
		
		// for each of the three waveguides	
		for(i=0; i<=2; ++i)
		{
			// find the fundamental mode, save it
			WMM_findfundmode(wg[i], 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[i], 100, 150, '0'+i, '0', 'C');
			// :-)
			m.fancymfile(fc, display[i], 100, 150, '0'+i, '0');
		}
	}
	return 0;
}