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

/*
 * 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.0;   // 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  30.0    // half length of the directional couplers, 
                      //    cores are considered to be decoupled for |z| > Z
#define Numz  30      // number of discretization steps in the numerical 
                      //    integration of the coupled mode equations

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

/* Parameter scans: range and stepsize */
#define Pmin  0.5     
#define Pmax  2.5
#define Pstep 0.05
// adjust the statement in line 248 !

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

/* 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;

	if(CalcBasisModes == 1)
	{
		wg = swgdef();
		WMM_findfundmode(wg, POL, SYM, 0.0, 0.0, par, 's', 't', ma);
		stmo = ma(0);
		stmo.write_def('s', 't');
		ma.clear();

		wg = rwgdef();
		WMM_findfundmode(wg, POL, SYM, 0.0, 0.0, par, 'r', 'i', ma);
		rimo = ma(0);
		rimo.write_def('r', 'i');
		ma.clear();
	}
	else 
	{
		stmo.read_def(POL, SYM, 's', 't');
		rimo.read_def(POL, SYM, 'r', 'i');
	}

	Dmatrix sigma(2,2);
	Dvector sbeta(2);
	Dmatrix samp(2,2);

	Dmatrix S(2,2);
	Dmatrix BK(2,2);

	double betaq = 0.5*(stmo.beta + rimo.beta);

	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);

	Cvector inA(2); inA(0) = CC1; inA(1) = CC0;
	Cvector ina(2); ina(0) = CC0; ina(1) = CC1;
	Cvector out;

	Cmatrix T(2,2);
	Cmatrix Tp(2,2);
	double h = WgpZ*2.0/((double)Numz);

for(double p=Pmin; p<=Pmax; p+=Pstep)
{
	// ... 
	Wgpb = p; 
	// ... 
	fprintf(stderr, "\n(( %g ))\n", p);

	T.unit(2);

	fprintf(stderr, "<< CMT: ");
	double z = -WgpZ;
	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();
	M0 = mult(add(mult(One, betaq), mult(mult(S, BK), -1.0)), CCI);  

	for(int zi=0; zi<=Numz-1; ++zi)
	{
		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);

		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);

		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);

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

	T = mult(T, expi(-betaq*(2.0*WgpZ)));
	fprintf(stderr, "T:\n");
	T.write(stderr);
	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", kappa);
	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", tau);

	fprintf(stderr, "|tau|^2+|kappa|^2: %g\n", tau+kappa);
	fprintf(stderr, "\n");

	apptoxyf("kappa", p, kappa);
	apptoxyf("tau", p, tau);
	apptoxyf("total", p, tau+kappa);
}
	ma.clear();

	return 0;
}