/* 
 * WMM-Application:
 * directional coupler with varying gap width, coupled mode theory
 */

/*
 * WMM
 * Wave-matching method for mode analysis of dielectric waveguides
 * Manfred Hammer 
 * (2002)
 */

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

/* coupler parameters (lengths in micrometers) */

#define Wgps   2.00   // straight waveguide: core width
#define Wgpv   0.14   // straight waveguide: core thickness
#define Wgpng  1.98   // straight waveguide: core refractive index
#define Wgpnb  1.45   // substrate/straight waveguide, background refr. index 

#define Wgpw   2.75   // ring waveguide: core width
#define Wgpt   1.0    // ring waveguide: core thickness
#define Wgpnr  1.596  // ring waveguide: core refractive index
#define Wgpna  1.0    // cover refractive index 

double Wgpb =  1.0;   // straight waveguide: burying depth 
double Wgpg =  0.5;   // straight waveguide: center displacement from ring rim 
double Wgpp =  0.0;   // position of ring waveguide; 0: closest approach 

#define Wgpl  1.55    // vacuum wavelength 
#define WgpR  100.0   // ring radius, outer rim

#define WgpZ  40.0    // half length of the directional couplers, 
                      //    cores are considered to be decoupled for |z| > Z
#define Numz  40      // number of discretization steps in the numerical 
                      //    integration of the coupled mode equations

Polarization POL=QTE; // QTE or QTM (semivectorial simulation)

#define CalcBasisModes 1 // 1: calculate the basis fields, save them to files
                         // 0: read them from previously stored files

/* display window, mode profile plots; x in [x0, x1], y in [y0, y1], 
   format: Disp(x0, y0, x1, y1), note the axis orientation! */
Rect Disp(-3.5, -4.9, 2.0, 3.9);

/* the straight waveguide */
Waveguide swgdef()
{
	Waveguide g(1, 1);

	g.hx(0) = 0.0;
	g.hx(1) = Wgpv;
	g.hy(0) = -Wgps/2.0;
	g.hy(1) =  Wgps/2.0;

	g.n(0,0) = Wgpnb;
	g.n(0,1) = Wgpnb;
	g.n(0,2) = Wgpnb;
	g.n(1,0) = Wgpnb;
	g.n(1,1) = Wgpng;
	g.n(1,2) = Wgpnb;
	g.n(2,0) = Wgpnb;
	g.n(2,1) = Wgpnb;
	g.n(2,2) = Wgpnb;

	g.lambda = Wgpl;

	return g;
}

/* the ring waveguide */
Waveguide rwgdef()
{
	Waveguide g(1, 1);

	g.hx(0) = 0.0;
	g.hx(1) = Wgpt;
	g.hy(0) = -Wgpw/2.0;
	g.hy(1) =  Wgpw/2.0;

	g.n(0,0) = Wgpnb;
	g.n(0,1) = Wgpnb;
	g.n(0,2) = Wgpnb;
	g.n(1,0) = Wgpna;
	g.n(1,1) = Wgpnr;
	g.n(1,2) = Wgpna;
	g.n(2,0) = Wgpna;
	g.n(2,1) = Wgpna;
	g.n(2,2) = Wgpna;

	g.lambda = Wgpl;

	return g;
}

/* the entire coupler structure */
#define Ltol 1.0e-8
Waveguide cpwgdef()
{
	if(Wgpnb <= Ltol)
	{
		fprintf(stderr, "mres.c: cpwgdef: Wgpb<=0 not implemented.\n");
		exit(1);
	}
	
	Dvector hytmp(4);
	hytmp(0) = -Wgpp-Wgpw; 
	hytmp(1) = -Wgpp; 
	hytmp(2) =  Wgpg-Wgps/2.0; 
	hytmp(3) =  Wgpg+Wgps/2.0; 
	hytmp.sort(0);
	Dvector hy(4);
	hy(0) = hytmp(0); 
	int ny=1;
	for(int j=1; j<=3; ++j)
	{
		if(fabs(hy(ny-1) - hytmp(j)) >= Ltol)
		{
			hy(ny) = hytmp(j); 
			++ny;
		}
	}
	--ny;
	Waveguide wg(3, ny);

	wg.hx(0) = -Wgpb-Wgpv;
	wg.hx(1) = -Wgpb;
	wg.hx(2) = 0.0;
	wg.hx(3) = Wgpt;
	
	for(int i=0; i<=ny; ++i) wg.hy(i) = hy(i);
	
	for(int m=0; m<=ny+1; ++m) 
	{
		double y;
		Rect r;
		r = wg.rectbounds(1, m);
		y = (r.y0+r.y1)/2.0;

		wg.n(0,m) = Wgpnb;
		if(y < Wgpg-Wgps/2.0 || y > Wgpg+Wgps/2.0) wg.n(1, m) = Wgpnb;
		else                                       wg.n(1, m) = Wgpng; 
		wg.n(2,m) = Wgpnb;
		if(y < -Wgpp-Wgpw || y > -Wgpp) wg.n(3, m) = Wgpna;
		else                            wg.n(3, m) = Wgpnr; 
		wg.n(4,m) = Wgpna;
	}

	wg.lambda = Wgpl;
	return wg;
}

/* mode solver parameters */
WMM_Parameters pardef()
{
	WMM_Parameters p;

	p.vform = HXHY;
	p.vnorm = NRMMH;
	p.ccomp = CCALL;

	p.ini_d_alpha   = 0.01;
	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;
}

/* coupled mode theory simulation */
int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray ma;
	WMM_Mode m, stmo, rimo;
	Waveguide wg;
	double bs, ba, lc;
	int n;

	// for field profile plots: select the interesting field component 
	Fcomp fcp = principalcomp(POL);

	if(CalcBasisModes == 1)
	{
		// define the straight waveguide
		wg = swgdef();
		// find its fundamental mode
		WMM_findfundmode(wg, POL, SYM, 0.0, 0.0, par, 's', 't', ma);
		// save it to a file
		stmo = ma(0);
		stmo.write_def('s', 't');
		ma.clear();

		// define the ring waveguide
		wg = rwgdef();
		// find its fundamental mode
		WMM_findfundmode(wg, POL, SYM, 0.0, 0.0, par, 'r', 'i', ma);
		// save it to a file
		rimo = ma(0);
		rimo.write_def('r', 'i');
		ma.clear();
	}
	else 
	{
		// read both basis modes from files 
		stmo.read_def(POL, SYM, 's', 't');
		rimo.read_def(POL, SYM, 'r', 'i');
	}

// look at the position of closest approach, z=0:

	// shift the basic modes to the positions of the two waveguides, 
	// compose an array of basis modes  
	m = stmo;
	m.translate(-Wgpb-Wgpv, Wgpg);
	ma.add(m);	
	m = rimo;
	m.translate(0.0, -(Wgpp+Wgpw/2.0));
	ma.add(m);	
	n = 2;

	// the coupler structure
	Wgpp = 0.0;
	wg = cpwgdef();

	// control output: basis waveguides and composite structures, geometry
	fprintf(stderr, "\n Straight waveguide (WG1):\n");
	ma(0).wg.write(stderr);
	fprintf(stderr, "\n Ring waveguide (WG2):\n");
	ma(1).wg.write(stderr);
	fprintf(stderr, "\n Composite coupler structure:\n");
	wg.write(stderr);

	// plots of the basis fields, shifted
	fprintf(stderr, "\nBasis mode profiles:\n"); 
	ma(0).mfile(fcp, ORG, Disp, 80, 50, 's', 't', 'S');
	ma(0).mfile(fcp, MOD, Disp, 120, 120, 's', 't', 'I');
	ma(1).mfile(fcp, ORG, Disp, 80, 50, 'r', 'i', 'S');
	ma(1).mfile(fcp, MOD, Disp, 120, 120, 'r', 'i', 'I');
	fprintf(stderr, "\n"); 

	// initiate coupled mode theory, 
	// compute supermode propagation constants and amplitude vectors
	Dmatrix sigma(n,n);
	Dvector sbeta(n);
	Dmatrix samp(n,n);

	Dmatrix S(n,n);
	Dmatrix BK(n,n);
	CMTmats(ma, wg, S, BK);
	
	// control output ...
	// fprintf(stderr, "[S]\n");
	// S.write(stderr);
	// fprintf(stderr, "[B+K]\n");
	// BK.write(stderr);
	
	// supermodes related to the z-invariant structure:  
	CMTsetup(ma, wg, sigma, sbeta, samp);

	bs = sbeta(0);
	ba = sbeta(1);
	// supermode propagation constants
	fprintf(stderr, "\nSupermodes in z=0:\n"); 
	fprintf(stderr, "beta_s: %g mum^(-1), vec(%g, %g)\n", 
		bs, samp(0,0), samp(1, 0)); 
	fprintf(stderr, "beta_a: %g mum^(-1), vec(%g, %g)\n", 
		ba, samp(0,1), samp(1, 1)); 
	// the coupling length
	lc = PI/fabs(bs-ba);
	fprintf(stderr, "   L_c: %g mum\n", lc); 
	fprintf(stderr, "\n"); 

	// compose the supermode profiles
	WMM_ModeArray smode;
	CMTsupermodes(ma, wg, sbeta, samp, smode);

	// plots of the supermode profiles
	fprintf(stderr, "\nSupermode profiles:\n"); 
	smode(0).mfile(fcp, ORG, Disp, 80, 50, 's', '0', 'S');
	smode(1).mfile(fcp, ORG, Disp, 80, 50, 'a', '0', 'S');
	smode(0).mfile(fcp, MOD, Disp, 120, 120, 's', '0', 'I');
	smode(1).mfile(fcp, MOD, Disp, 120, 120, 'a', '0', 'I');
	ma.clear();

// ring/straight waveguide coupling, longitudinally varying structure

	// memory for output curves
	Dvector pAB(Numz);
	Dvector pAb(Numz);
	Dvector paB(Numz);
	Dvector pab(Numz);
	Dvector len(Numz);
	len = Dvector(Numz);

	// an average propagation constant 
	double betaq = 0.5*(stmo.beta + rimo.beta);

	// matrices involved in the Runge-Kutta integration
	Cmatrix M0(2,2);
	Cmatrix Mp(2,2);
	Cmatrix M1(2,2);
	Cmatrix N1(2,2);
	Cmatrix N2(2,2);
	Cmatrix N3(2,2);
	Cmatrix N4(2,2);
	Cmatrix One;
	One.unit(2);

	// initial mode amplitudes 
	// input in coupler port A
	Cvector inA(2); inA(0) = CC1; inA(1) = CC0;
	// input in coupler port a
	Cvector ina(2); ina(0) = CC0; ina(1) = CC1;
	// local mode amplitudes  
	Cvector out;

	// stepsize for the numerical integration
	double h = WgpZ*2.0/((double)Numz);
	// coupler transfer matrix, contribution from an integration step
	Cmatrix Tp(2,2);
	// cummulative coupler transfer matrix
	Cmatrix  T(2,2);

	fprintf(stderr, "\n<< CMT: ");
	// start with T=1 at z=-Z
	T.unit(2);
	double z = -WgpZ;
	Wgpp = WgpR-sqrt(WgpR*WgpR-z*z);
	ma.clear();

	// for properly shifted cores ...
	m = stmo; m.translate(-Wgpb-Wgpv, Wgpg); ma.add(m);	
	m = rimo; m.translate(0.0, -(Wgpp+Wgpw/2.0)); ma.add(m);	
	wg = cpwgdef();
	// ... evaluate the CMT matrices ...
	CMTmats(ma, wg, S, BK);
	S.inverse();
	// ... for amplitudes with the explicit average phase evolution 
	M0 = mult(add(mult(One, betaq), mult(mult(S, BK), -1.0)), CCI);  

	for(int zi=0; zi<=Numz-1; ++zi)
	{
		// half a stepsize advanced
		z += h/2.0;
		Wgpp = WgpR-sqrt(WgpR*WgpR-z*z);
		ma.clear();
		m = stmo; m.translate(-Wgpb-Wgpv, Wgpg); ma.add(m);	
		m = rimo; m.translate(0.0, -(Wgpp+Wgpw/2.0)); ma.add(m);	
		wg = cpwgdef();
		CMTmats(ma, wg, S, BK);
		S.inverse();
		Mp = mult(add(mult(One, betaq), mult(mult(S, BK), -1.0)), CCI);

		// another half stepsize advanced
		z += h/2.0;
		Wgpp = WgpR-sqrt(WgpR*WgpR-z*z);
		ma.clear();
		m = stmo; m.translate(-Wgpb-Wgpv, Wgpg); ma.add(m);	
		m = rimo; m.translate(0.0, -(Wgpp+Wgpw/2.0)); ma.add(m);	
		wg = cpwgdef();
		CMTmats(ma, wg, S, BK);
		S.inverse();
		M1 = mult(add(mult(One, betaq), mult(mult(S, BK), -1.0)), CCI);
		
		// apply the Runge-Kutta recipe	
		N1 = mult(M0, h);
		N2 = mult(mult(Mp, add(One, mult(N1, 0.5))), h);
		N3 = mult(mult(Mp, add(One, mult(N2, 0.5))), h);
		N4 = mult(mult(M1, add(One, N3)), h);
		Tp = One;
		Tp = add(Tp, mult(N1, 1.0/6.0));
		Tp = add(Tp, mult(N2, 1.0/3.0));
		Tp = add(Tp, mult(N3, 1.0/3.0));
		Tp = add(Tp, mult(N4, 1.0/6.0));
		
		T = mult(Tp, T);

		
		// evaluate the local amplitudes: remove the average phase 
		Tp = mult(T, expi(-betaq*(z+WgpZ)));

		// store the local amplitudes
		len(zi) = z;
		out = mult(Tp, inA); 
		pAB(zi) = out(0).sqabs(); 
		pAb(zi) = out(1).sqabs();
		out = mult(Tp, ina); 
		paB(zi) = out(0).sqabs(); 
		pab(zi) = out(1).sqabs();

		M0 = M1;
		fprintf(stderr, ".");
	}
	fprintf(stderr, " Ok. >>\n");

	// the scattering matrix of the entire coupler
	T = mult(T, expi(-betaq*(2.0*WgpZ)));
	fprintf(stderr, "\nCoupler scattering matrix T:\n\n");
	T.write(stderr);
	// coupling constant: both offdiagonal elements should be equal
	fprintf(stderr, "|kappa|^2: %g\n", T(1, 0).sqabs());
	fprintf(stderr, "           %g\n", T(0, 1).sqabs());
	double kappa = 0.5*(T(0, 1).sqabs()+T(1, 0).sqabs());
	fprintf(stderr, "  average: %g\n\n", kappa);
	// transfer constant: both diagonal elements should have equal
	//                    absolute values
	double tau = T(0, 0).sqabs();
	double rho = T(1, 1).sqabs();
	fprintf(stderr, "|tau|^2: %g\n", tau);
	fprintf(stderr, "|rho|^2: %g\n", rho);	
	tau = 0.5*(tau+rho);
	fprintf(stderr, "average: %g\n\n", tau);

	// power conservation: |tau|^2+|kappa|^2 = 1 ?
	fprintf(stderr, "|tau|^2+|kappa|^2: %g\n\n", tau+kappa);
	
	// save amplitude evolution to files
	fprintf(stderr, "Amplitude evolution, excitation in A:\n");
	toxyf("pAB", len, pAB);	
	toxyf("pAb", len, pAb);	
	toxyf("pAB+pAb", len, add(pAB, pAb));	
	fprintf(stderr, "Amplitude evolution, excitation in a:\n");
	toxyf("paB", len, paB);	
	toxyf("pab", len, pab);	
	toxyf("paB+pab", len, add(paB, pab));	
	fprintf(stderr, "\nOk.\n\n");

	ma.clear();

	return 0;
}