/* 
 * METRIC-Application:
 * Rectangular slab waveguide facet
 */

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

#define Pol    TE    // light polarization
#define FPnb   1.0   // background refractive index
#define FPng   2.0   // core refractive index
#define FPw    0.4   // core width /mum
#define Wavel  1.55  // vacuum wavelength /mum

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

#define DWx    2.5   // display window, distances from the facet center 
#define DWz    3.0   //

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

/* the facet */
Circuit facetdef()
{
	Circuit fac(1, 1);
	fac.hx(0) = -FPw/2.0;
	fac.hx(1) =  FPw/2.0;
	fac.hz(0) = -FPw;
	fac.hz(1) =  0.0;
	fac.n(0, 0) = FPnb; fac.n(0, 1) = FPnb; fac.n(0, 2) = FPnb;
	fac.n(1, 0) = FPng; fac.n(1, 1) = FPng; fac.n(1, 2) = FPnb;
	fac.n(2, 0) = FPnb; fac.n(2, 1) = FPnb; fac.n(2, 2) = FPnb;
	fac.lambda = Wavel;
	return fac;
}

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

/* facet simulation */
int main()
{
	fprintf(stderr, "\nWaveguide facet:\n");	
	fprintf(stderr, "----------------\n");	

	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", FPnb);
	fprintf(stderr, "n_g    = %g\n", FPng);
	fprintf(stderr, "w      = %g mum\n\n", FPw);

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

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

	fprintf(stderr, "\nReflected fundamental mode amplitude:\n"); 
	Complex a; 
	double p;
	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);
	
	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);
	fprintf(stderr, "Total scattered guided power:\n");
	fprintf(stderr, "P_back    = %g\n\n", pb);

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

	qfld.adjustphase(0.0, -FPw/2.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; 
}