/* 
 * WMM-Application:
 * polarization splitter based on radiatively coupled waveguides 
 * coupled mode theory
 */

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

/* coupler parameters */
#define Wgpns 1.95  // substrate refractive index
#define Wgpnf 2.2   // film refractive index 
#define Wgpnc 1.0   // cover: air
#define Wgpl  1.3   // vacuum wavelength 
#define Wgpw  1.21  // outer rib width
#define Wgph  0.7   // rib height
#define WgpW  4.96  // width of the central rib
#define Wgpg  0.11  // width of the gaps
#define WgpL  850.0 // optimum device length


/* one single port waveguides */
Waveguide swgdef()
{
	Waveguide g(1, 1);
	g.hx(0) = 0.0;
	g.hx(1) = 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) = 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;
}

/* the central rib */
Waveguide cwgdef()
{
	Waveguide g(1, 1);
	g.hx(0) = 0.0;
	g.hx(1) = 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) = 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;
}

/* the three waveguide coupler structure */
Waveguide rcwgdef()
{
	Waveguide g(1, 5);
	g.hx(0) = 0.0;
	g.hx(1) = Wgph;
	g.hy(0) = -Wgpw-Wgpg-WgpW/2.0;
	g.hy(1) =      -Wgpg-WgpW/2.0;
	g.hy(2) =           -WgpW/2.0;
	g.hy(3) =            WgpW/2.0;
	g.hy(4) =       Wgpg+WgpW/2.0;
	g.hy(5) =  Wgpw+Wgpg+WgpW/2.0;
	g.n(0,0) = Wgpns;
	g.n(0,1) = Wgpns;
	g.n(0,2) = Wgpns;
	g.n(0,3) = Wgpns;
	g.n(0,4) = Wgpns;
	g.n(0,5) = Wgpns;
	g.n(0,6) = Wgpns;
	g.n(1,0) = Wgpnc;
	g.n(1,1) = Wgpnf;
	g.n(1,2) = Wgpnc;
	g.n(1,3) = Wgpnf;
	g.n(1,4) = Wgpnc;
	g.n(1,5) = Wgpnf;
	g.n(1,6) = Wgpnc;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = Wgpnc;
	g.n(2,2) = Wgpnc;
	g.n(2,3) = Wgpnc;
	g.n(2,4) = Wgpnc;
	g.n(2,5) = Wgpnc;
	g.n(2,6) = Wgpnc;
	g.lambda = Wgpl;
	return g;
}

/* two port waveguides at a distance */
Waveguide pwgdef()
{
	Waveguide g(1, 3);
	g.hx(0) = 0.0;
	g.hx(1) = Wgph;
	g.hy(0) = -Wgpw-Wgpg-WgpW/2.0;
	g.hy(1) =      -Wgpg-WgpW/2.0;
	g.hy(2) =       Wgpg+WgpW/2.0;
	g.hy(3) =  Wgpw+Wgpg+WgpW/2.0;
	g.n(0,0) = Wgpns;
	g.n(0,1) = Wgpns;
	g.n(0,2) = Wgpns;
	g.n(0,3) = Wgpns;
	g.n(0,4) = Wgpns;
	g.n(1,0) = Wgpnc;
	g.n(1,1) = Wgpnf;
	g.n(1,2) = Wgpnc;
	g.n(1,3) = Wgpnf;
	g.n(1,4) = Wgpnc;
	g.n(2,0) = Wgpnc;
	g.n(2,1) = Wgpnc;
	g.n(2,2) = Wgpnc;
	g.n(2,3) = Wgpnc;
	g.n(2,4) = Wgpnc;
	g.lambda = Wgpl;
	return g;
}

/* analysis 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   = 20;
	p.ini_alpha_max = 2.0;
	p.ini_steps = 50;
	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 = 2.5;
	p.btol   = 1.0e-7;
	p.mshift = 1.0e-8;
	return p;
}

/* polarization splitting ratios */
double er0(double p0te, double p0tm)
{
	return 10.0*log10(p0te/p0tm);
}
double er1(double p1te, double p1tm)
{
	return 10.0*log10(p1tm/p1te);
}

int main()
{
	WMM_Parameters par = pardef();
	Polarization pol;
	int m, svp;
	WMM_ModeArray psi;
	WMM_ModeArray cmodes;
	WMM_Mode mode;
	WMM_Mode som;
	Waveguide wg;
	int ii, o0, o1;
	char c;
	int zi;
	double z;
	double neffmax = 0.0;
// space for power curves, P_1,2^TE,TM(z)
	#define Lmin    0.0 
	#define Lmax 1500.0 
	#define NumL 1500 
	Dmatrix p0(2, NumL+1);  
	Dmatrix p1(2, NumL+1);  
// space for power curves, P_1,2^TE,TM(delta q)
	#define Numq 200 
	Dmatrix wp0(2, Numq+1);
	Dmatrix wp1(2, Numq+1);
	Dmatrix Wp0(2, Numq+1);
	Dmatrix Wp1(2, Numq+1);
	Dmatrix gp0(2, Numq+1);
	Dmatrix gp1(2, Numq+1);
	Dmatrix hp0(2, Numq+1);
	Dmatrix hp1(2, Numq+1);
	Dmatrix lp0(2, Numq+1);
	Dmatrix lp1(2, Numq+1);

// for TE and TM polarization
	for(svp=0; svp<=1; ++svp)
	{
		if(svp==0) pol = QTE;
		else       pol = QTM;
		if(pol == QTE) { neffmax = 2.16; }
		if(pol == QTM) { neffmax = 2.11; }
		c = poltochar2(pol);
		fprintf(stderr, "(((T%c)))\n", c);

// if not available from former computation: compute basic mode set,
// all modes of the center waveguide, save to files
		cmodes.clear();
		wg = cwgdef();
		WMM_modeanalysis(wg,pol,SYM,Wgpns,neffmax,par,'c','-',psi);
		cmodes.merge(psi);
		psi.clear();
		WMM_modeanalysis(wg,pol,ASY,Wgpns,neffmax,par,'c','-',psi);
		cmodes.merge(psi);
		psi.clear();
		cmodes.write_def('t', c);
		cmodes.clear();
// the mode of the outer waveguide, save to file
		wg = swgdef();
		WMM_findfundmode(wg,pol,SYM,Wgpns,neffmax,par,'c','-',psi);
		psi(0).write_def('o', '0');
		psi.clear();

// second run should start here 
// load basic modes
		cmodes.read_def('t', c);
		som.read_def(pol, SYM, 'o', '0');
// shift som to port positions, add to coupled mode set
		mode = som;
		mode.translate(0.0, -(WgpW/2.0+Wgpg+Wgpw/2.0));
		cmodes.add(mode);
		mode = som;
		mode.translate(0.0, +(WgpW/2.0+Wgpg+Wgpw/2.0));
		cmodes.add(mode);
// m coupled modes
		m = cmodes.num;
// power input by mode with index ii
		ii = m-2;
// power detection in modes with indices o0, o1
		o0 = m-2;
		o1 = m-1;
// define the entire structure 
		wg = rcwgdef();
// initialize coupled mode theory
		Dmatrix sigma(m,m);
		Dvector sbeta(m);
		Dmatrix samp(m,m);
		CMTsetup(cmodes, wg, sigma, sbeta, samp);
// compute power transfer from channels ii to o0,o1 for varying device length
		for(zi = 0; zi<=NumL; ++zi)
		{
			z = Lmin+(Lmax-Lmin)/NumL*zi; 
			p0(svp, zi) = CMTpower(ii, o0, z, sigma, sbeta, samp);
			p1(svp, zi) = CMTpower(ii, o1, z, sigma, sbeta, samp);
		}
// save to files
		char pnam[20] = "p0_t_";
		Dvector zv(NumL+1);
		Dvector pv(NumL+1);
		for(zi = 0; zi<=NumL; ++zi) zv(zi) = Lmin+(Lmax-Lmin)/NumL*zi; 
		pnam[2] = '0'; pnam[4] = c; 
		for(zi = 0; zi<=NumL; ++zi) pv(zi) = p0(svp, zi); 
		toxyf(pnam, zv, pv);
		pnam[2] = '1'; pnam[4] = c; 
		for(zi = 0; zi<=NumL; ++zi) pv(zi) = p1(svp, zi); 
		toxyf(pnam, zv, pv);
// display window 
		Rect disp(-Wgph/2.0, -WgpW/2.0-Wgpg-Wgpw-Wgpw/2,
		          Wgph*1.5,  WgpW/2.0+Wgpg+Wgpw+Wgpw/2);
// compose CMT supermodes, save modes and field profiles to disk
		WMM_ModeArray smodes;
		CMTsupermodes(cmodes, wg, sbeta, samp, smodes);
		for(m=0; m<=smodes.num-1; ++m) smodes(m).write_def('S', '0'+m);
		Fcomp fc;
		fc = principalcomp(pol); 
		for(m=0; m<=smodes.num-1; ++m)
		      smodes(m).mfile(fc, SQR, disp, 70, 150, 'S', '0'+m, 'C');
// make a propagation plot
		// the outlets
		WMM_ModeArray pmodes;
		pmodes.clear();
		pmodes.add(cmodes(o0));
		pmodes.add(cmodes(o1));
		Waveguide pwg;
		pwg = pwgdef();
		// one mode is excited 
		Cvector iamp(2);
		iamp(0) = CC1;
		iamp(1) = CC0;
		// no perturbation
		Cvector pert(cmodes.num);
		pert.init(CC0); 
		// write to matlab file
		smodes.ioprop(iamp, pmodes, pwg, pmodes, pwg, pert, 
		              Wgph/2.0, WgpL, 
		              disp.y0, disp.y1, 80,
		              -200.0, WgpL+200.0, 200, 
		              't', c);
		smodes.clear();

// for tolerancing: store gradients of the supermode propagation constants
// with respect to geometry parameters
		int numm;
		WMM_Mode mode;
		numm = cmodes.num;
		cmodes.clear();
		Dvector gradW(numm);
		Dvector gradw(numm);
		Dvector gradg(numm);
		Dvector gradh(numm);
		Dvector gradl(numm);
		for(m=0; m<=numm-1; ++m)
		{
			mode.read_def(pol, NOS, 'S', '0'+m);
			gradW(m) = 0.0;
			gradW(m) += mode.vergeovar(1, 3);
			gradW(m) += mode.vergeovar(1, 4);
			gradW(m) += mode.vergeovar(1, 5);
			fprintf(stderr, "(%d) grad W: %g\n", m, gradW(m));
			gradw(m) = 0.0;
			gradw(m) += mode.vergeovar(1, 5);
			gradw(m) -= mode.vergeovar(1, 0);
			fprintf(stderr, "(%d) grad w: %g\n", m, gradw(m));
			gradg(m) = 0.0;
			gradg(m) += mode.vergeovar(1, 4);
			gradg(m) += mode.vergeovar(1, 5);
			gradg(m) -= mode.vergeovar(1, 1);
			gradg(m) -= mode.vergeovar(1, 0);
			fprintf(stderr, "(%d) grad g: %g\n", m, gradg(m));
			gradh(m) = 0.0;
			gradh(m) += mode.horgeovar(1, 1);
			gradh(m) += mode.horgeovar(1, 3);
			gradh(m) += mode.horgeovar(1, 5);
			fprintf(stderr, "(%d) grad h: %g\n", m, gradh(m));
			gradl(m) = -(mode.beta + Wgpw*gradw(m)
					       + WgpW*gradW(m) 
					       + Wgpg*gradg(m) 
					       + Wgph*gradh(m))/Wgpl;
		}

// for tolerancing: store output power for original mode profiles, but with
// perturbed propagation constants, save to files
		int qi;
		double dq;
		Dvector qsbeta(numm);
		for(qi=0; qi<=Numq; ++qi)
		{
			dq = -0.1+0.2/Numq*qi;
			for(m=0; m<=numm-1; ++m) 
				qsbeta(m) = sbeta(m) + gradw(m)*dq;
			wp0(svp,qi) = CMTpower(ii,o0,WgpL,sigma,qsbeta,samp);
			wp1(svp,qi) = CMTpower(ii,o1,WgpL,sigma,qsbeta,samp);
			for(m=0; m<=numm-1; ++m) 
				qsbeta(m) = sbeta(m) + gradW(m)*dq;
			Wp0(svp,qi) = CMTpower(ii,o0,WgpL,sigma,qsbeta,samp);
			Wp1(svp,qi) = CMTpower(ii,o1,WgpL,sigma,qsbeta,samp);
			for(m=0; m<=numm-1; ++m) 
				qsbeta(m) = sbeta(m) + gradg(m)*dq;
			gp0(svp,qi) = CMTpower(ii,o0,WgpL,sigma,qsbeta,samp);
			gp1(svp,qi) = CMTpower(ii,o1,WgpL,sigma,qsbeta,samp);
			for(m=0; m<=numm-1; ++m) 
				qsbeta(m) = sbeta(m) + gradh(m)*dq;
			hp0(svp,qi) = CMTpower(ii,o0,WgpL,sigma,qsbeta,samp);
			hp1(svp,qi) = CMTpower(ii,o1,WgpL,sigma,qsbeta,samp);
			for(m=0; m<=numm-1; ++m) 
				qsbeta(m) = sbeta(m) + gradl(m)*dq;
			lp0(svp,qi) = CMTpower(ii,o0,WgpL,sigma,qsbeta,samp);
			lp1(svp,qi) = CMTpower(ii,o1,WgpL,sigma,qsbeta,samp);
		}
		char qnam[20] = "_p0_t_";
		Dvector qv(Numq+1);
		Dvector rv(Numq+1);
		for(qi = 0; qi<=Numq; ++qi) qv(qi) = -0.1+0.2/Numq*qi;
		qnam[0] = 'w'; qnam[3] = '0'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = wp0(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'w'; qnam[3] = '1'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = wp1(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'W'; qnam[3] = '0'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = Wp0(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'W'; qnam[3] = '1'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = Wp1(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'g'; qnam[3] = '0'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = gp0(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'g'; qnam[3] = '1'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = gp1(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'h'; qnam[3] = '0'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = hp0(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'h'; qnam[3] = '1'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = hp1(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'l'; qnam[3] = '0'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = lp0(svp, qi);
		toxyf(qnam, qv, rv); 
		qnam[0] = 'l'; qnam[3] = '1'; qnam[5] = c; 
		for(qi = 0; qi<=Numq; ++qi) rv(qi) = lp1(svp, qi);
		toxyf(qnam, qv, rv); 
	}

// characterization of polarization splitting needs data for TE and TM
// evaluate polarization splitting ratios versus length, save 
	char enam[20] = "er_";
	Dvector zv(NumL+1);
	Dvector ev(NumL+1);
	for(zi = 0; zi<=NumL; ++zi) zv(zi) = Lmin+(Lmax-Lmin)/NumL*zi; 
	enam[2] = '0'; 
	for(zi = 0; zi<=NumL; ++zi) ev(zi) = er0(p0(0, zi), p0(1, zi));
	toxyf(enam, zv, ev);
	enam[2] = '1'; 
	for(zi = 0; zi<=NumL; ++zi) ev(zi) = er1(p1(0, zi), p1(1, zi));
	toxyf(enam, zv, ev);

// for tolerancing:
// evaluate polarization splitting ratios versus geometry paramters w, W, g, h,
// write to files
	int qi;
	char qnam[20] = "_er_";
	Dvector qv(Numq+1);
	Dvector rv(Numq+1);
	for(qi = 0; qi<=Numq; ++qi) qv(qi) = -0.1+0.2/Numq*qi;
	qnam[0] = 'w'; qnam[3] = '0';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er0(wp0(0, qi), wp0(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'w'; qnam[3] = '1';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er1(wp1(0, qi), wp1(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'W'; qnam[3] = '0';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er0(Wp0(0, qi), Wp0(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'W'; qnam[3] = '1';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er1(Wp1(0, qi), Wp1(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'g'; qnam[3] = '0';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er0(gp0(0, qi), gp0(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'g'; qnam[3] = '1';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er1(gp1(0, qi), gp1(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'h'; qnam[3] = '0';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er0(hp0(0, qi), hp0(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'h'; qnam[3] = '1';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er1(hp1(0, qi), hp1(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'l'; qnam[3] = '0';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er0(lp0(0, qi), lp0(1, qi));
	toxyf(qnam, qv, rv); 
	qnam[0] = 'l'; qnam[3] = '1';
	for(qi = 0; qi<=Numq; ++qi) rv(qi) = er1(lp1(0, qi), lp1(1, qi));
	toxyf(qnam, qv, rv); 

	return 0;
}