/* 
 * WMM-Application:
 * coupling between two rib waveguides 
 */

#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 Wgpg  0.5  // gap width
#define Wgpl  1.3  // vacuum wavelength 

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

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

	g.n(0,0) = Wgpns;
	g.n(0,1) = Wgpns;
	g.n(0,2) = Wgpns;
	g.n(0,3) = Wgpns;
	g.n(0,4) = Wgpns;
	g.n(1,0) = Wgpnf;
	g.n(1,1) = Wgpnf;
	g.n(1,2) = Wgpnf;
	g.n(1,3) = Wgpnf;
	g.n(1,4) = Wgpnf;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = Wgpnf;
	g.n(2,2) = Wgpnc;
	g.n(2,3) = Wgpnf;
	g.n(2,4) = Wgpnc;
	g.n(3,0) = Wgpnc;
	g.n(3,1) = Wgpnc;
	g.n(3,2) = Wgpnc;
	g.n(3,3) = Wgpnc;
	g.n(3,4) = 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.1;
	p.ini_N_alpha   = 10;
	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 QTE modes,
   save modes, 
   write mode profile data, 
   write interference data */
int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	WMM_Mode m;
	int i;
	double lc;
		
	// define the waveguide
	Waveguide wg = wgdef();

	// find its two fundamental mode of opposite symmetry, save them
	WMM_findfundmode(wg, VEC, ASY, 0.0, 0.0, par, '-', '-', ma);
	m = ma(0);
	WMM_findfundmode(wg, VEC, SYM, 0.0, 0.0, par, '-', '-', ma);
	ma.add(m);
	ma.write_def('-', '-');

	// second run should start here
//	ma.read_def('-', '-');

//	Mlop_Print=YES;

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

	// make plots of the mode profiles,
	for(i=0; i<=ma.num-1; ++i)
		ma(i).mfile(EY, ORG, display, 95, 125, '0', '0', 'S'); 

	// no perturbation
	Cvector p(ma.num);
	p.init(CC0);

	// equal initial mode amplitudes 
	Cvector a(ma.num);
	a.init(CC1);
	
	// interference animation, frames for two coupling lengths 
	lc = PI/(ma(0).beta - ma(1).beta);	

	ma.movie(a, p, EY, MOD, display, 75, 105, 'I', 20, 0.0, 2.0*lc);
	ma.movie(a, p, EY, MOD, display, 75, 105, 'F', 20, 0.0, 2.0*lc);

	return 0;
}