/* 
 * METRIC-Application:
 * T-junction of slab waveguides
 */

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

#define Pol    TE    // light polarization
#define TPnb   1.0   // background refractive index
#define TPng   3.4   // core refractive index
#define TPw    0.22  // core width /mum
#define Wavel  1.55  // vacuum wavelength /mum

#define CWx    7.0   // computational window, 
#define CWz    7.0   // distances from the coner center 
#define NumMx  100   // number of spectral terms in the field discretization
#define NumMz  100   //

#define DWx    4.0   // display window, distances from the corner center 
#define DWz    4.0   //

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

/* the junction configuration */
Circuit cornerdef()
{
	Circuit cr(1, 1);
	cr.hx(0) = -TPw/2.0;
	cr.hx(1) =  TPw/2.0;
	cr.hz(0) = -TPw/2.0;
	cr.hz(1) =  TPw/2.0;
	cr.n(2, 0) = TPnb; cr.n(2, 1) = TPng; cr.n(2, 2) = TPnb;
	cr.n(1, 0) = TPng; cr.n(1, 1) = TPng; cr.n(1, 2) = TPnb;
	cr.n(0, 0) = TPnb; cr.n(0, 1) = TPng; cr.n(0, 2) = TPnb;
	cr.lambda = Wavel;
	return cr;
}

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

/* simulation of the T-junction */
int main()
{
	fprintf(stderr, "\nWaveguide T-junction:\n");	
	fprintf(stderr,   "---------------------\n");	

	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", TPnb);
	fprintf(stderr, "n_g    = %g\n", TPng);
	fprintf(stderr, "w      = %g mum\n\n", TPw);

	// the structure under investigation 
	Circuit cr = cornerdef();
	// vertical computational window
	Interval vcw(-CWx, CWx);
	// horizontal computational window
	Interval hcw(-CWz, CWz);
	// QUEP setup 
	QuepField qfld(cr, Pol, vcw, NumMx, hcw, NumMz);

	// input: the fundamental mode in the left channel 
	qfld.input(LEFT, CC1, 0);
	// calculate the light propagation along the junction
	qfld.quepsim(); 

	Complex a; 
	double p;
	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n"); 
	a = qfld.Aout(LEFT, 0);
	p = qfld.Pout(LEFT, 0);
	fprintf(stderr, "R = %g + i %g,  |R|^2 = %g\n", a.re, a.im, p);
	a = qfld.Aout(TOP, 0);
	p = qfld.Pout(TOP, 0);
	fprintf(stderr, "U = %g + i %g,  |U|^2 = %g\n", a.re, a.im, p);
	a = qfld.Aout(BOTTOM, 0);
	p = qfld.Pout(BOTTOM, 0);
	fprintf(stderr, "D = %g + i %g,  |D|^2 = %g\n", a.re, a.im, p);
	
	double pb, pf, pd, pu;
	pb = qfld.Pout(LEFT);
	pf = qfld.Pout(RIGHT);
	pd = qfld.Pout(BOTTOM);
	pu = qfld.Pout(TOP);
	
	fprintf(stderr, "Total scattered power:\n");
	fprintf(stderr, "P_back    = %g\n", pb);
	fprintf(stderr, "P_forw    = %g\n", pf);
	fprintf(stderr, "P_down    = %g\n", pd);
	fprintf(stderr, "P_up      = %g\n", pu);
	fprintf(stderr, "Sum P_out = %g\n", pb+pf+pd+pu);
	
	pb = qfld.Pgout(LEFT);
	pu = qfld.Pgout(TOP);
	pd = qfld.Pgout(BOTTOM);
	fprintf(stderr, "Total scattered guided power:\n");
	fprintf(stderr, "P_back    = %g\n", pb);
	fprintf(stderr, "P_up      = %g\n", pu);
	fprintf(stderr, "P_down    = %g\n\n", pd);

	char pc = polchr(Pol);
	Fcomp fc = principalcomp(Pol);
	Mlop_Print = NO;
	Mlop_Colour = YES;

	qfld.adjustphase(0.0, -TPw/16.0);
	qfld.viewer(-DWx, DWx, -DWz, DWz, 150, 150, 't', pc); 
	qfld.plot(fc, REP, -DWx, DWx, -DWz, DWz, 150, 150, 't', pc); 
	qfld.plot(fc, IMP, -DWx, DWx, -DWz, DWz, 150, 150, 't', pc); 
	qfld.plot(fc, MOD, -DWx, DWx, -DWz, DWz, 150, 150, 't', pc);
	qfld.phasemap(fc, -DWx, DWx, -DWz, DWz, 150, 150, 't', pc);
	qfld.movie(-DWx, DWx, -DWz, DWz, 150, 150, 30, 't', pc);             
	qfld.fplot(fc, -DWx, DWx, -DWz, DWz, 150, 150, 't', pc); 
	qfld.fmovie(-DWx, DWx, -DWz, DWz, 150, 150, 30, 't', pc);             

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