/* 
 * METRIC-Application:
 * two adjacent parallel waveguides, 
 * conventional coupled mode theory model 
 */

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

#define Pol TM
#define MPnb       1.45     /* background refractive index */ 
#define MPng       3.40     /* guiding regions: refractive index */ 
#define Wavel      1.55     /* vacuum wavelength / mum */
double GPw  = 0.20;         /* width of the waveguides / mum */
double GPg  = 0.30;         /* gap between the waveguides / mum */

/* Scan over gap widths: parameter interval */
#define MinG 0.1           
#define MaxG 0.9           
#define NumG 100

/* display window, animation of field propagation */
#define DWx0  -1.5
#define DWx1   1.5
#define DWz0   0.0
#define DWz1  15.0

/* ---------------------------------------------------------- */ 

/* a single waveguide */
Waveguide iowgdef()
{
	Waveguide g(1);
	g.hx(0) = -GPw/2.0;
	g.hx(1) =  GPw/2.0;
	g.n(0) = MPnb;
	g.n(1) = MPng;
	g.n(2) = MPnb;
	g.lambda = Wavel;
	return g;
}

/* the two parallel waveguides */
Waveguide cpwgdef()
{
        Waveguide g(3);
        g.hx(0) = -(GPg/2.0+GPw);
        g.hx(1) = -(GPg/2.0);
        g.hx(2) =   GPg/2.0;
        g.hx(3) =   GPg/2.0+GPw;
        g.n(0) = MPnb;
        g.n(1) = MPng;
        g.n(2) = MPnb;
        g.n(3) = MPng;
        g.n(4) = MPnb;
        g.lambda = Wavel;
        return g;
}                     

/* ---------------------------------------------------------- */ 

/* something went wrong: stop the program */
void cmtstop(const char *e)
{
	fprintf(stderr, "CMT: ERROR: %s.\n", e);
	exit(1);
}

/* ---------------------------------------------------------- */ 

/* coupled waveguide simulation */
int main()
{
	int j;
	char pc = polchr(Pol);
	Fcomp fcp = principalcomp(Pol);

	ModeArray ma, basis;
	Mode m_m, m_p;

	// the input/output waveguides and the coupling segment
	Waveguide wg, cp;

	fprintf(stderr, "\nTwo parallel waveguides:\n");	
	fprintf(stderr, "----------------------\n\n");	

	fprintf(stderr, "Device parameters:\n");	
	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", MPnb);
	fprintf(stderr, "n_g    = %g\n", MPng);
	fprintf(stderr, "w   = %g mum\n", GPw);
	fprintf(stderr, "g   = %g mum\n", GPg);
	fprintf(stderr, "\n");	

	// display window for mode plots
        Interval display(-GPg-2.0*GPw, GPg+2.0*GPw);       

/* modes of the port waveguides */
	wg = iowgdef();
	fprintf(stderr, "Port modes:\n");

	modeanalysis(wg, Pol, ma);
	if(ma.num != 1) cmtstop("Input guide is not single mode");
	m_m = ma(0);
	fprintf(stderr, "T%c0:  beta = %g /mum,  n_eff = %g\n", 
		pc, m_m.beta, m_m.neff);
	m_m.translate(-GPg/2.0-GPw/2.0);
	m_m.plot(fcp, ORG, display, 500, 'p', 'm', 'V');
	m_p = ma(0);
	m_p.translate(GPg/2.0+GPw/2.0);
	m_p.plot(fcp, ORG, display, 500, 'p', 'p', 'V');

/* the CMT mode basis */
	basis.clear();
	basis.add(m_m);
	basis.add(m_p);

/* compute supermode propagation constants and amplitude vectors */
	Dmatrix sigma;
	Dvector smpc;
	Dmatrix smav;
	cp = cpwgdef();
	CMTsetup(basis, cp, sigma, smpc, smav);

	fprintf(stderr,"Sigma:\n");
	sigma.write(stderr);
	fprintf(stderr,"Smpc:\n");
	smpc.write(stderr);
	fprintf(stderr,"Smav:\n");
	smav.write(stderr);

	fprintf(stderr, "\nCMT-supermodes:\n");
	fprintf(stderr, "    T%c0:  beta = %g /mum,  n_eff = %g\n", 
			pc, smpc(0), smpc(0)/2.0/PI*Wavel);
	fprintf(stderr, "    T%c1:  beta = %g /mum,  n_eff = %g\n", 
			pc, smpc(1), smpc(1)/2.0/PI*Wavel);
	fprintf(stderr, "    CMT: L_c = %g mum\n", PI/(smpc(0)-smpc(1)));


/* The exact solution: */
	fprintf(stderr, "\nExact coupling:\n");
	cp = cpwgdef();

	fprintf(stderr, "\nExact supermodes:\n");
	modeanalysis(cp, Pol, ma);
	for(j=0; j<=ma.num-1; ++j) 
	{
		fprintf(stderr, "    T%c%d:  beta = %g /mum,  n_eff = %g\n", 
			pc, j, ma(j).beta, ma(j).neff);
		char ct = 'V'; if(j==1) ct = 'L';
		ma(j).plot(fcp, ORG, display, 500, 'r', dig1(j), ct);
	}
	fprintf(stderr, "    Exact: L_c = %g mum\n", 
			PI/(ma(0).beta-ma(1).beta));
	fprintf(stderr, "\n");

/* Field animations */
	ModeArray sma;
	Cvector amp_f, amp_b;

	// first waveguide
	sma.clear();
	sma.add(basis(0));
	amp_f = Cvector(sma.num); amp_f.init(CC0);
	amp_b = Cvector(sma.num); amp_b.init(CC0);
	amp_f(0) = CC1;
	sma.plot(amp_f,amp_b,fcp,REP,DWx0,DWx1,DWz0,DWz1,80,200,'s','0');
	sma.movie(amp_f,amp_b,DWx0,DWx1,DWz0,DWz1,80,200,30,'s','0');

	// second waveguide
	sma.clear();
	sma.add(basis(1));
	amp_f = Cvector(sma.num); amp_f.init(CC0);
	amp_b = Cvector(sma.num); amp_b.init(CC0);
	amp_f(0) = CC1;
	sma.plot(amp_f,amp_b,fcp,REP,DWx0,DWx1,DWz0,DWz1,80,200,'s','1');
	sma.movie(amp_f,amp_b,DWx0,DWx1,DWz0,DWz1,80,200,30,'s','1');

	// the coupled waveguides 
	amp_f = Cvector(ma.num); amp_f.init(CC0);
	amp_b = Cvector(ma.num); amp_b.init(CC0);
	amp_f(0) = CC1;
	amp_f(1) = (-1)*CC1;
	ma.plot(amp_f,amp_b,fcp,REP,DWx0,DWx1,DWz0,DWz1,80,200,'c','p');
	ma.movie(amp_f,amp_b,DWx0,DWx1,DWz0,DWz1,80,200,30,'c','p');

	fprintf(stderr, "\n");
	fprintf(stderr, "Coupling length versus gap width:\n");
	Dvector gval(NumG+1);
	Dvector cmtlc(NumG+1);
	Dvector exalc(NumG+1);

	for(int gi=0; gi<=NumG; ++gi)
	{
		GPg = MinG+((double)gi)/((double)NumG)*(MaxG-MinG);
		gval(gi) = GPg;
		
		wg = iowgdef();
		modeanalysis(wg, Pol, ma, 1);
		m_m = ma(0);
		m_m.translate(-GPg/2.0-GPw/2.0);
		m_p = ma(0);
		m_p.translate(GPg/2.0+GPw/2.0);
		basis.clear();
		basis.add(m_m);
		basis.add(m_p);
		cp = cpwgdef();
		CMTsetup(basis, cp, sigma, smpc, smav);
		cmtlc(gi) = PI/(smpc(0)-smpc(1));

		modeanalysis(cp, Pol, ma, 1);
		exalc(gi) = PI/(ma(0).beta-ma(1).beta);
	}
	fprintf(stderr, "\nLc_cmt   vs. g >>>");
	toxyf("clc", gval, cmtlc);
	fprintf(stderr, "Lc_exact vs. g >>>");
	toxyf("elc", gval, exalc);

	fprintf(stderr, "\nOkay. \n");
	return 0; 
}