/* 
 * WMM-Application:
 * MMI coupler, 1x2/2x1 power splitter/combiner 
 */

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

/* waveguide parameters */
#define Wgpns 1.95  // substrate refractive index
#define Wgpnf 2.2   // film refractive index
#define Wgpnc 1.0   // cover: air
#define Wgpw  1.3   // input waveguide width 
#define WgpW  4.45  // width of the multimode segment 
#define Wgph  0.7   // strip height
#define Wgpl  1.3   // vacuum wavelength 
#define WgpS  1.05  // off-center position of port waveguides
#define WgpL  32.3  // length of the MMI segment 

/* waveguide definition: a single strip of width w */
Waveguide wgdef(double w)
{
	Waveguide g(1, 1);

	g.hx(0) = 0.0;
	g.hx(1) = Wgph;
	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;
}

/* two output ports */
Waveguide wg2def()
{
	Waveguide g(1, 3);

	g.hx(0) = 0.0;
	g.hx(1) = Wgph;
	g.hy(0) = -Wgpw/2.0-WgpS;
	g.hy(1) =  Wgpw/2.0-WgpS;
	g.hy(2) = -Wgpw/2.0+WgpS;
	g.hy(3) =  Wgpw/2.0+WgpS;

	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 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   = 15;
	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 = 3.0;

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

	return p;
}

/* calculate QTE modes, save modes, 
   write mode profile data, 
   write interference data */
int main()
{
	WMM_Parameters par = pardef();
	WMM_ModeArray mmis;
	WMM_ModeArray mmia;
	WMM_ModeArray mmi;
	WMM_Mode m;
	WMM_Mode t;
	Waveguide wg;
	int i;

	// the input/output waveguide
	wg = wgdef(Wgpw);
	// calculate the input/output mode
	WMM_findfundmode(wg, QTE, SYM, 1.95, 2.11, par, '-', '-', mmi);
	// ... and save it
	m = mmi(0);
	m.write_def('-', '-');
	mmi.clear();
	
	// the multimode segment 
	wg = wgdef(WgpW);
	// find all its symmetrical modes
	WMM_modeanalysis(wg, QTE, SYM, 1.95, 2.11, par, '-', '-', mmis);
	// and save to files 
	mmis.write_def('s', '-');
	// find all its antriysymmetrical modes
	WMM_modeanalysis(wg, QTE, ASY, 1.95, 2.11, par, '-', '-', mmia);
	// and save to files 
	mmia.write_def('a', '-');
	
	// second run should start here ...
//	m.read_def(QTE, SYM, '-', '-');
//	mmis.read_def('s', '-');
//	mmia.read_def('a', '-');

	// display window
	Rect display(-Wgph, -(WgpW/2.0+1.0), Wgph*1.5, WgpW/2.0+1.0);

//	Mlop_Print=YES;
//	Mlop_Colour=NO;
	
	// write mode profiles, surface and contour plots
	m.mfile(EY, ORG, display, 75, 155, '0', 'i', 'S'); 
	m.mfile(EY, SQR, display, 75, 155, '0', 'i', 'C'); 
	for(i=0; i<=mmis.num-1; ++i)
	{
		mmis(i).mfile(EY,ORG,display,75,155,dig10(i),dig1(i),'S'); 
		mmis(i).mfile(EY,SQR,display,75,155,dig10(i),dig1(i),'C'); 
	}
	for(i=0; i<=mmia.num-1; ++i)
	{
		mmia(i).mfile(EY,ORG,display,75,155,dig10(i),dig1(i),'S'); 
		mmia(i).mfile(EY,SQR,display,75,155,dig10(i),dig1(i),'C'); 
	}

	// no perturbation
	Cvector p(mmis.num+mmia.num);
	p.init(CC0);

	// first we need only symmetrical modes	
	mmi = mmis;

	// initial mode amplitudes: centered input waveguide 
	Cvector a(mmi.num);
	for(i=0; i<=mmi.num-1; ++i)
	{
		t = mmi(i);
		a(i) = Complex(lscalprod(t, m), 0.0);
	}
	// --> the initial intensity distribution
	mmi.mfile(a, p, 0.0, SZ, ORG, display, 75, 155, '0', 's', 'C');

	// propagation along one focusing length
	mmi.prop(a, p, Wgph/2.0, display.y0, display.y1, 115,
                                 0.0, WgpL, 250, 'p', '0');

	// interference animation 
//	mmi.movie(a, p, EY, MOD, display, 75, 155, 'I', 32, 0.0, WgpL);
//	mmi.movie(a, p, SZ, ORG, display, 75, 155, 'F', 32, 0.0, WgpL);

	// to determine a proper device length:
	// power transfer versus length of the MMI segment
        // centered input port -> centered output port
	mmi.writeiopower(m, m, p, 1000, 0.0, 100.0, '0', '0');

	// to determine the off-center port position:
	// transverse intensity contours at half the focusing length
	mmi.mfile(a, p, WgpL/2.0, SZ, ORG, display, 75, 155, '1', 's', 'C');


	// interference: top views, port waveguides included

        // one complete focus length, single centered input/output ports 
	mmi.ioprop(m, m, p, Wgph/2.0, WgpL, display.y0, display.y1, 115, 
                   -2.0, WgpL+2.0, 250, '0', '0');

	// power splitter  
	WMM_ModeArray omode;
	omode.add(m);
	omode.add(m);
	omode(0).translate(0.0, -WgpS);
	omode(1).translate(0.0,  WgpS);
	Waveguide owg=wg2def();
	mmi.ioprop(m, omode, owg, p, Wgph/2.0, WgpL/2.0, 
                   display.y0, display.y1, 115, 
                   -2.0, WgpL/2.0+2.0, 250, 's', '0');
	omode.clear();

	// power combiner 
	// equal imput amplitudes, constructive interference  
	WMM_ModeArray imode;
	imode.add(m);
	imode.add(m);
	imode(0).translate(0.0, -WgpS);
	imode(1).translate(0.0,  WgpS);
	Waveguide iwg=wg2def();
	Cvector iamp(2);
	iamp(0) =  CC1;
	iamp(1) =  CC1;
	omode.add(m);
	mmi.ioprop(iamp, imode, iwg,  omode, m.wg, p, Wgph/2.0, WgpL/2.0, 
                   display.y0, display.y1, 115, 
                   -2.0, WgpL/2.0+2.0, 250, 'c', '0');
	// 180^o phase difference of the input modes, destructive interference  
	// now the antisymmetrical modes are required as well
	mmi.merge(mmia);
	mmi.sort();
	iamp(0) =  CC1;
	iamp(1) = -CC1;
	mmi.ioprop(iamp, imode, iwg,  omode, m.wg, p, Wgph/2.0, WgpL/2.0, 
                   display.y0, display.y1, 115, 
                   -2.0, WgpL/2.0+2.0, 250, 'd', '0');
	// only one of the input channels of the power combiner is excited: 
	// about half of the input power passes the device
	iamp(0) =  CC1;
	iamp(1) =  CC0;
	mmi.ioprop(iamp, imode, iwg,  omode, m.wg, p, Wgph/2.0, WgpL/2.0, 
                   display.y0, display.y1, 115, 
                   -2.0, WgpL/2.0+2.0, 250, 'a', '0');

	return 0;
}

