/*
 * WMM-Application:
 * nonreciprocal TE - phaseshifter,
 * rib waveguide, polar magnetooptic configuration,
 * compensation wall / domain lattice  
 */

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


/* rib waveguide parameters */
#define Wgpt       0.5 // film thickness, after etching (in mum)
#define Wgpw       1.0 // width of the waveguide (in mum)
#define Wgph       0.5 // etching depth (in mum)
#define Wgpnf      2.2 // film refractive index 
#define Wgpns      1.9 // substrate refractive index
#define Wgpnc      1.0 // cover: air
#define Wgpl       1.3 // vacuum wavelength (in mum)
#define WgpFrot 3147.0 // specific Faraday-rotation in ^o/cm 

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

	g.hx(0) = 0.0;
	g.hx(1) = Wgpt;
	g.hx(2) = 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) = 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   = 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 shifts */
int main()
{
	WMM_Parameters par = pardef();
	int nfm;
	WMM_ModeArray ma;
	WMM_Mode mode;
	int i;
	Complex db;
	Perturbation p;
	double vz;

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

	// calculate its fundamental semivectorial TE mode
	nfm = WMM_findfundmode(wg, QTE, SYM, 0.0, 0.0, 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(-(Wgph+Wgpt),-2.0*Wgpw,1.3*(Wgph+Wgpt),2.0*Wgpw);
		mode.mfile(EY, MOD, display, 75, 115, '0', '0', 'C'); 

	// calculate the nonreciprocal phase shift for a
	// polar domain lattice, lattice period = rib width
		db = CC0;
		// two domains inside the rib 
		p = magopt_polar(-WgpFrot, Wgpnf, Wgpl, 
				 Rect(0.0, 0.0, Wgpt+Wgph, Wgpw/2.0));
		db = db + mode.phaseshift(p);
		p = magopt_polar( WgpFrot, Wgpnf, Wgpl, 
				 Rect(0.0, -Wgpw/2.0, Wgpt+Wgph, 0.0));
		db = db + mode.phaseshift(p);
		// several domains in the lateral film
		vz = -1.0;
		for(i=1; i<=20; ++i)
		{
			p = magopt_polar(-vz*WgpFrot, Wgpnf, Wgpl, 
		            Rect(0.0, i*Wgpw/2.0, Wgpt, (i+1)*Wgpw/2.0));
			db = db + mode.phaseshift(p);
			p = magopt_polar(vz*WgpFrot, Wgpnf, Wgpl, 
			    Rect(0.0, -(i+1)*Wgpw/2.0, Wgpt, -i*Wgpw/2.0));
			db = db + mode.phaseshift(p);
			vz *= -1.0;
		}
		// the difference between the propagation constants of
		// forward and backward propagating modes 
		fprintf(stderr, "dl-NRPS:  %5g cm^(-1)\n", 2.0*db.re*10000.0);

	// calculate the nonreciprocal phase shift for a
	// compensation wall at the rib center, separating two polar domains  
		// perturbations on the rib area ...
		db = CC0; 
		p = magopt_polar(-WgpFrot, Wgpnf, Wgpl, 
				 Rect(0.0, 0.0, Wgpt+Wgph, Wgpw/2.0));
		db = db + mode.phaseshift(p);
		p = magopt_polar( WgpFrot, Wgpnf, Wgpl, 
				 Rect(0.0, -Wgpw/2.0, Wgpt+Wgph, 0.0));
		db = db + mode.phaseshift(p);
		// ... and outside the rib
		p = magopt_polar(-WgpFrot, Wgpnf, Wgpl, 
				 Rect(0.0, Wgpw/2.0, Wgpt, 10.0));
		db = db + mode.phaseshift(p);
		p = magopt_polar( WgpFrot, Wgpnf, Wgpl, 
				 Rect(0.0, -10.0, Wgpt, -Wgpw/2.0));
		db = db + mode.phaseshift(p);
		// the difference between the propagation constants of
		// forward and backward propagating modes 
		fprintf(stderr, "cw-NRPS:  %g cm^(-1)\n", 2.0*db.re*10000.0);
	}
	ma.clear();

	return 0;
}