/*
 * WMM-Application:
 * nonreciprocal TM - phaseshifter,
 * rib waveguide, equatorial magnetooptic configuration
 * double layer magnetooptic film
 */

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


/* rib waveguide parameters */
#define Wgpt         0.48 // total lateral film thickness (in mum)
#define Wgpb         0.20 // thickness of the bottom magnetooptic layer
#define Wgph         0.04 // rib height, etching depth (in mum)
#define Wgpw         1.5  // rib width  (in mum)
#define Wgpns        1.95 // substrate refractive index
#define Wgpnp        2.27 // refractive index of the bottom layer
#define WgpFrotp   350.0  // bottom layer: specific Faraday-rotation in ^o/cm 
#define Wgpnn        2.33 // refractive index of the top layer
#define WgpFrotn -1450.0  // top layer: specific Faraday-rotation in ^o/cm 
#define Wgpnc        1.0  // cover: air
#define Wgpl         1.3  // vacuum wavelength (in mum)

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

	g.hx(0) = 0.0;
	g.hx(1) = Wgpb;
	g.hx(2) = Wgpt;
	g.hx(3) = Wgpt+Wgph;
	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) = Wgpnp;
	g.n(1,1) = Wgpnp;
	g.n(1,2) = Wgpnp;
	g.n(2,0) = Wgpnn;
	g.n(2,1) = Wgpnn;
	g.n(2,2) = Wgpnn;
	g.n(3,0) = Wgpnc;
	g.n(3,1) = Wgpnn;
	g.n(3,2) = Wgpnc;
	g.n(4,0) = Wgpnc;
	g.n(4,1) = Wgpnc;
	g.n(4,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   = 10;
	p.ini_alpha_max = 2.0;

	p.ini_steps = 20;
	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-8;
	return p;
}

/* calculate fundamental mode and nonreciprocal phase shift */
int main()
{
	WMM_Parameters par = pardef();
	int nfm;
	WMM_ModeArray ma;
	WMM_Mode mode;
	Complex db;
	Perturbation p;

	// define the waveguide 
	Waveguide wg = wgdef();

	// calculate its fundamental semivectorial TM mode
	nfm = WMM_findfundmode(wg, QTM, SYM, Wgpns, Wgpnn, par, '-', '-', ma);
	if(nfm >= 1)
	{
		mode = ma(0);
		// save the mode
		mode.write_def('0','0');
		// make a contour plot of the mode profile
		Rect display(-2.0*Wgpt,-3.0*Wgpw,1.3*(Wgph+Wgpt),3.0*Wgpw);
		mode.mfile(HY, MOD, display, 75, 115, '0', '0', 'C'); 

		// calculate the nonreciprocal phase shift
		db = CC0;
		// contribution of the bottom layer 
		p = magopt_equat(WgpFrotp, Wgpnp, Wgpl, 
				 Rect(0.0, -10.0, Wgpb, 10.0));
		db = db + mode.phaseshift(p);
		// contribution of the top layer 
		p = magopt_equat(WgpFrotn, Wgpnn, Wgpl, 
				 Rect(Wgpb, -10.0, Wgpt, 10.0));
		db = db + mode.phaseshift(p);
		p = magopt_equat(WgpFrotn, Wgpnn, Wgpl, 
				 Rect(Wgpt, -Wgpw/2.0, Wgpt+Wgph, Wgpw/2.0));
		db = db + mode.phaseshift(p);
		// the difference between the propagation constants of
		// forward and backward propagating modes 
		fprintf(stderr, "NRPS:  %5g cm^(-1)\n", 2.0*db.re*10000.0);
	}
	ma.clear();

	return 0;
}