/* 
 * WMM-Application:
 * Unidirectional polarization converter: 
 * double layer magnetooptic raised strip waveguide, 
 * TE/TM mode conversion, coupled mode theory model  
 */

/*
 * WMM
 * Wave-matching method for mode analysis of dielectric waveguides
 * Manfred Lohmeyer 
 * University of Osnabrueck, Department of Physics
 * Barbarastrasse 7, D-49069 Osnabrueck, Germany
 * (1999)
 */

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

/* rib waveguide parameters */
#define Wgpnf      2.302 // film refractive index
#define Wgpns      1.95  // substrate refractive index
#define Wgpnc      1.0   // cover: air
#define Wgpl       1.3   // vacuum wavelength [mum]
#define WgpFrot 3007.42  // specific Faraday rotation of the two layers (+/-)   
                         // [degrees/cm]


/* rasied strip waveguide, heigth h, width w */
Waveguide wgdef(double h, double w)
{
	Waveguide g(1, 1);

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

	g.n(0,0) = Wgpns;
	g.n(0,1) = Wgpns;
	g.n(0,2) = Wgpns;
	g.n(1,0) = Wgpnc;
	g.n(1,1) = Wgpnf;
	g.n(1,2) = Wgpnc;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = Wgpnc;
	g.n(2,2) = Wgpnc;

	g.lambda = Wgpl;

	return g;
}

/* Mode amplitudes according to coupled mode theory:
   two coupled modes, Hermitian coupling matrix, 
   b0, b1: propagation constants,
   kappa:  coupling coefficient, 
   c0, c1: initial amplitudes, 	 
   z:      propagation distance
   j:      0,1; index of the returned amplitude */
Complex cmtcoeff(int j, double b0, double b1, Complex kappa, 
		Complex c0, Complex c1, double z)
{
	double db, rho;
	Complex m0, m1, a;
	db = b0-b1;
	rho = sqrt(db*db/4.0+kappa.sqabs());

	switch(j)
	{
	case 0:
		m0 = Complex(cos(rho*z), -db/2.0/rho*sin(rho*z));
		m1 = kappa*(-sin(rho*z)/rho);
	        m1 = m1*CCI;
		break;
	case 1:
		m0 = kappa.conj()*(-sin(rho*z)/rho);
	        m0 = m0*CCI;
		m1 = Complex(cos(rho*z), db/2.0/rho*sin(rho*z));
	}
	a = m0*c0+m1*c1;	
	a = a*expi(-(b0+b1)/2.0*z);
	return a;
}

/* 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   = 20;
	p.fin_alpha_max = 2.0;

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

	return p;
}

/* unidirectional polarization converters: */
int main()
{

// four different geometries: find the fundamental modes, save them  

	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	double w, h;
	int i;

	// strip height and width
	double hdata[4] = {0.6341, 0.8260, 1.0118, 1.1975};
	double wdata[4] = {0.7976, 0.9951, 1.1914, 1.3865};

	for(i=0; i<=3; i+=1)
	{
		w = wdata[i];
		h = hdata[i];
		Waveguide wg = wgdef(h, w);

		WMM_findfundmode(wg, VEC, SYM, 0.0, 0.0, par, '-', '-', ma);
		ma(0).write_def('o', 'a'+i);
		ma.clear();

		WMM_findfundmode(wg, VEC, ASY, 0.0, 0.0, par, '-', '-', ma);
		ma(0).write_def('o', 'a'+i);
		ma.clear();
	}

// magnetooptic perturbation: 
// coupled mode theory for two fundamental modes

	WMM_Mode tem;
	WMM_Mode tmm;
	double pte, ptm;
	double ktete, ktetm, ktmtm;
	double kte, ktm, kap;
	double t, theta, xi;

	// bottom layer thicknesses 
	double wgptdata[4] = {0.1807, 0.2643, 0.3491, 0.4311};
	// magnetization orientation [degrees]
	double wgpthdata[4] = {62.8, 65.0, 66.2, 67.9};
			
	// for the four geometries	
	for(i=0; i<=3; i+=1)
	{
		fprintf(stderr, "\n\n");
		
	// load the unperturbed modes
		tem.read_def(VEC, SYM, 'o', 'a'+i);
		tmm.read_def(VEC, ASY, 'o', 'a'+i);

	// show the actual geometry
		w = tem.wg.hy(1) - tem.wg.hy(0);
		h = tem.wg.hx(1) - tem.wg.hx(0);
		t = wgptdata[i];
		theta = wgpthdata[i];
		fprintf(stderr, "[%c]  h = %g  w = %g  t = %g  theta = %g\n", 
			         'a'+i, h, w, t, theta);

	// modes are assumed to be normalized
		pte = 1.0;
		ptm = 1.0;

	// two regions with opposite specific Faraday rotation
		Rect rI(t, -w/2.0, h, w/2.0);     // top layer
		Rect rII(0.0, -w/2.0, t, w/2.0);  // bottom layer 

	// compute the entries of the coupling matrix
		xi = Wgpnf*Wgpl/PI*WgpFrot/180.0*PI/10000.0;
		ktete = twomrecint(tem, EX, tem, EZ, rI)
		      - twomrecint(tem, EX, tem, EZ, rII);
		ktete *= xi/val_invomep0(tem.wg.lambda)
			 /2.0/pte; 
		ktmtm = twomrecint(tmm, EX, tmm, EZ, rI)
		      - twomrecint(tmm, EX, tmm, EZ, rII);
		ktmtm *= xi/val_invomep0(tem.wg.lambda)
			 /2.0/ptm; 
		ktetm = twomrecint(tem, EX, tmm, EY, rI)
		      - twomrecint(tem, EY, tmm, EX, rI)
		      - twomrecint(tem, EX, tmm, EY, rII)
		      + twomrecint(tem, EY, tmm, EX, rII);
		ktetm *= xi/val_invomep0(tem.wg.lambda)
			 /sqrt(pte)/sqrt(ptm)/4.0; 
			
		t = theta/180.0*PI;
		kte = ktete*sin(t);
		ktm = ktmtm*sin(t);
		kap = ktetm*cos(t);

	// phase mismatch and coupling coefficient for forward and
	// backward propagation, [cm^(-1)], conversion length
		fprintf(stderr, "F: b_te+k_te-b_tm-k_tm = %g\n",
			(tem.beta+kte-tmm.beta-ktm)*10000.0);
		fprintf(stderr, "F: kappa               = i %g\n",
			kap*10000.0);
		fprintf(stderr, "B: b_te-k_te-b_tm+k_tm = %g\n",
			(tem.beta-kte-tmm.beta+ktm)*10000.0);
		fprintf(stderr, "F: -kappa              = i %g\n",
			-kap*10000.0);
		fprintf(stderr, "L                      = %g\n",
			PI/2.0/kap);	

	// plots of the mode power conversion, forward and backward direction,
        // TE and TM polarization
		char ftename[20] = "fe_.xyf";
		char ftmname[20] = "fm_.xyf";
		char btename[20] = "be_.xyf";
		char btmname[20] = "bm_.xyf";
		FILE *fedat;
		FILE *fmdat;
		FILE *bedat;
		FILE *bmdat;
		ftename[2] = 'a'+i;
		ftmname[2] = 'a'+i;
		btename[2] = 'a'+i;
		btmname[2] = 'a'+i;
		fedat = fopen(ftename, "w+");
		fmdat = fopen(ftmname, "w+");
		bedat = fopen(btename, "w+");
		bmdat = fopen(btmname, "w+");

		int zi;
		double z;
		Complex c;
		double b0f, b1f, b0b, b1b; 
		b0f = tem.beta+kte;	
		b0b = tem.beta-kte;	
		b1f = tmm.beta+ktm;	
		b1b = tmm.beta-ktm;	
		for(zi = 0; zi<=500; ++zi)
		{
			z = zi*(PI/2.0/fabs(kap))/500;

			c = cmtcoeff(0,b0f,b1f,Complex(0.0,kap),CC1,CC0,z);
			fprintf(fedat, "%g  %g\n", z/1000.0, c.sqabs());

			c = cmtcoeff(1,b0f,b1f,Complex(0.0,kap),CC1,CC0,z);
			fprintf(fmdat, "%g  %g\n", z/1000.0, c.sqabs());

			c = cmtcoeff(0,b0b,b1b,Complex(0.0,-kap),CC1,CC0,z);
			fprintf(bedat, "%g  %g\n", z/1000.0, c.sqabs());

			c = cmtcoeff(1,b0b,b1b,Complex(0.0,-kap),CC1,CC0,z);
			fprintf(bmdat, "%g  %g\n", z/1000.0, c.sqabs());
		}
		fclose(fedat);
		fclose(fmdat);
		fclose(bedat);
		fclose(bmdat);
	}
	return 0;
}
