/* 
 * METRIC-Application:
 * A resonant waveguide junction, enabled by a square cavity
 */

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

#define Pol    TE     // light polarization
#define JPnb   1.00   // background refractive index
#define JPng   3.40   // waveguide cores: refractive index
#define JPw    0.100  // width of the port waveguides /mum
#define JPgh   0.355  // gap between lower port and cavity /mum 
                      // = gap between lower port and end of right port
#define JPgv   0.385  // gap between right port and cavity /mum
#define JPW    1.786  // width of the square cavity /mum
#define Wavel  1.55   // vacuum wavelength /mum

#define CWx    5.0    // computational window 
#define CWz    5.0    // (symmetrical in x, z, distance from the cavity center)
#define NumMx  120    // number of spectral terms in the field discretizations
#define NumMz  120    //

#define DWx    2.1    // display window 
#define DWz    2.1    // (symmetrical in x, z, distance from the cavity center)

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

/* the junction configuration */
Circuit sqrjuncdef()
{
	Circuit sqr;
	sqr = Circuit(3, 3);

	sqr.hx(0) = -JPW/2.0-JPgh-JPw;
	sqr.hx(1) = -JPW/2.0-JPgh;
	sqr.hx(2) = -JPW/2.0;
	sqr.hx(3) =  JPW/2.0;

	sqr.hz(0) = -JPW/2.0;
	sqr.hz(1) =  JPW/2.0;
	sqr.hz(2) =  JPW/2.0+JPgv;
	sqr.hz(3) =  JPW/2.0+JPgv+JPw;

	sqr.n(0, 0) = JPnb;
	sqr.n(0, 1) = JPnb;
	sqr.n(0, 2) = JPnb;
	sqr.n(0, 3) = JPnb;
	sqr.n(0, 4) = JPnb;

	sqr.n(1, 0) = JPng;
	sqr.n(1, 1) = JPng;
	sqr.n(1, 2) = JPng;
	sqr.n(1, 3) = JPng;
	sqr.n(1, 4) = JPng;

	sqr.n(2, 0) = JPnb;
	sqr.n(2, 1) = JPnb;
	sqr.n(2, 2) = JPnb;
	sqr.n(2, 3) = JPnb;
	sqr.n(2, 4) = JPnb;

	sqr.n(3, 0) = JPnb;
	sqr.n(3, 1) = JPng;
	sqr.n(3, 2) = JPnb;
	sqr.n(3, 3) = JPng;
	sqr.n(3, 4) = JPnb;

	sqr.n(4, 0) = JPnb;
	sqr.n(4, 1) = JPnb;
	sqr.n(4, 2) = JPnb;
	sqr.n(4, 3) = JPng;
	sqr.n(4, 4) = JPnb;

	sqr.lambda = Wavel;
	return sqr;
}

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

/* resonator simulation */
int main()
{
	fprintf(stderr, "\nSquare resonator + waveguide junction:\n");	
	fprintf(stderr, "--------------------------------------\n");	
	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", JPnb);
	fprintf(stderr, "n_g    = %g\n", JPng);
	fprintf(stderr, "w      = %g mum\n", JPw);
	fprintf(stderr, "g_h    = %g mum\n", JPgh);
	fprintf(stderr, "g_v    = %g mum\n", JPgv);
	fprintf(stderr, "W      = %g mum\n\n", JPW);

	Circuit sqr = sqrjuncdef();
	// vertical computational window
	Interval cwx(-CWx, CWx);
	// horizontal computational window
	Interval cwz(-CWz, CWz);
	// QUEP setup + field discretization 
	QuepField qfld(sqr, Pol, cwx, NumMx, cwz, NumMz);

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

	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n"); 
	Complex a;
	a = qfld.Aout(LEFT, 0);
	fprintf(stderr, "B = %g + i %g,  |B|^2 = %g\n", a.re, a.im, a.sqabs());
	a = qfld.Aout(RIGHT, 0);
	fprintf(stderr, "F = %g + i %g,  |F|^2 = %g\n", a.re, a.im, a.sqabs());
	a = qfld.Aout(TOP, 0);
	fprintf(stderr, "U = %g + i %g,  |U|^2 = %g\n", a.re, a.im, a.sqabs());
	
	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);
	pf = qfld.Pgout(RIGHT);
	pu = qfld.Pgout(TOP);
	fprintf(stderr, "Total scattered guided power:\n");
	fprintf(stderr, "P_back      = %g\n", pb);
	fprintf(stderr, "P_forw      = %g\n", pf);
	fprintf(stderr, "P_up        = %g\n", pu);
	fprintf(stderr, "Sum P_out,g = %g\n\n", pb+pf+pu);

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

	qfld.adjustphase();
	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; 
}