/* 
 * WMM-Application:
 * two parallel strips, coupled mode theory
 */

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

/* waveguide parameters */
#define Wgpw  1.2  // strip width
#define Wgph  0.7  // strip height
#define Wgpg  0.25 // gap width
#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 

/* a single waveguide */
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 entire coupler structure */
Waveguide cpwgdef()
{
	Waveguide g(1, 3);

	g.hx(0) = 0.0;
	g.hx(1) = Wgph;

	g.hy(0) = -Wgpw-Wgpg/2.0;
	g.hy(1) =      -Wgpg/2.0;
	g.hy(2) =       Wgpg/2.0;
	g.hy(3) =  Wgpw+Wgpg/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;
}

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

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

	// define one waveguide
	wg = swgdef();
	// find its fundamental mode
	WMM_findfundmode(wg, QTM, SYM, 0.0, 0.0, par, 'o', '-', ma);
	// save it to a file
	som = ma(0);
	som.write_def('o', '0');
	ma.clear();
	// read it again 
	// som.read_def(pol, SYM, 'o', '0');

	// shift it to the positions of the two waveguides, 
	// compose an array of coupled modes 
	m = som;
	m.translate(0.0, -(Wgpg/2.0+Wgpw/2.0));
	ma.add(m);	
	m = som;
	m.translate(0.0, +(Wgpg/2.0+Wgpw/2.0));
	ma.add(m);	
	n = 2;

	// the coupler structure
	wg = cpwgdef();

	// init coupled mode theory, 
	// compute supermode propagation constants and amplitude vectors
	Dmatrix sigma(n,n);
	Dvector sbeta(n);
	Dmatrix samp(n,n);
	CMTsetup(ma, wg, sigma, sbeta, samp);
	bs = sbeta(0);
	ba = sbeta(1);
	fprintf(stderr, "beta_s: %.15g mum^(-1)\n", bs); 
	fprintf(stderr, "beta_a: %.15g mum^(-1)\n", ba); 
	// the coupling length
	lc = PI/fabs(bs-ba);
	fprintf(stderr, "   L_c: %.15g mum\n", lc); 

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

	// surface plots of the mode profiles
	Rect disp(-Wgph, -Wgpg/2.0-Wgpw-Wgpw, 1.5*Wgph, Wgpg/2.0+Wgpw+Wgpw);
	smode(0).mfile(HY, ORG, disp, 95, 125, 's', '0', 'S');
	smode(1).mfile(HY, ORG, disp, 95, 125, 'a', '0', 'S');

// supermode interference:
	// equal initial amplitudes
	Cvector a(2);
	a(0) = CC1;
	a(1) = CC1;
	// no perturbation
	Cvector p(2);
	p(0) = CC0;
	p(1) = CC0;
	// intensity plot over five coupling lengths
	smode.prop(a, p, Wgph/3.0, disp.y0, disp.y1, 90, 
	                        0.0, 5.0*lc, 500, '0', '0');

	ma.clear();
	return 0;
}