/* 
 * METRIC-Application:
 * A rectangular microresonator, two parallel port waveguides 
 */

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

#define Pol   TE      // light polarization
#define RPnb  1.0     // background refractive index 
#define RPng  3.2     // guiding regions: refractive index 
#define RPw   0.2     // port waveguide width /mum 
#define RPg   0.29    // gap between waveguides and resonator /mum 
#define RPW   1.54    // cavity width /mum
#define RPL   1.54    // cavity length /mum
#define Wavel 1.532   // vacuum wavelength /mum

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

#define DWx  2.6      // vertical display window, dist. from the cavity center
#define DWz0 -2.0     // horizontal display window
#define DWz1 (RPL+2.0)

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

/* the resonator configuration */
Circuit resdef()
{
	Circuit s;
	s = Circuit(5, 1);
	s.hx(0) = -RPW/2.0-RPg-RPw;
	s.hx(1) = -RPW/2.0-RPg;
	s.hx(2) = -RPW/2.0;
	s.hx(3) =  RPW/2.0;
	s.hx(4) =  RPW/2.0+RPg;
	s.hx(5) =  RPW/2.0+RPg+RPw;
	s.hz(0) = 0.0;
	s.hz(1) = RPL;
	s.n(0, 0) = RPnb; s.n(0, 1) = RPnb; s.n(0, 2) = RPnb;
	s.n(1, 0) = RPng; s.n(1, 1) = RPng; s.n(1, 2) = RPng;
	s.n(2, 0) = RPnb; s.n(2, 1) = RPnb; s.n(2, 2) = RPnb;
	s.n(3, 0) = RPnb; s.n(3, 1) = RPng; s.n(3, 2) = RPnb;
	s.n(4, 0) = RPnb; s.n(4, 1) = RPnb; s.n(4, 2) = RPnb;
	s.n(5, 0) = RPng; s.n(5, 1) = RPng; s.n(5, 2) = RPng;
	s.n(6, 0) = RPnb; s.n(6, 1) = RPnb; s.n(6, 2) = RPnb;
	s.lambda = Wavel;
	return s;
}

/* a port waveguide */
Waveguide iowgdef()
{
        Waveguide g(1);
        g.hx(0) = -RPw/2.0;
        g.hx(1) =  RPw/2.0;
        g.n(0) = RPnb;
        g.n(1) = RPng;
        g.n(2) = RPnb;
        g.lambda = Wavel;
        return g;
}

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

/* resonator simulation */
int main()
{
	fprintf(stderr, "\nRectangular microresonator:\n");	
	fprintf(stderr, "---------------------------\n");	
	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", RPnb);
	fprintf(stderr, "n_g    = %g\n", RPng);
	fprintf(stderr, "w      = %g mum\n", RPw);
	fprintf(stderr, "g      = %g mum\n", RPg);
	fprintf(stderr, "W      = %g mum\n", RPW);
	fprintf(stderr, "L      = %g mum\n\n", RPW);

	// modes of the port waveguides 
	ModeArray ma;
	Mode m_m, m_p;	
	Waveguide iowg = iowgdef();
	fprintf(stderr, "Port modes:\n");
	modeanalysis(iowg, Pol, ma);
	m_m = ma(0);
	m_m.translate(-RPW/2.0-RPg-RPw/2.0);
	m_p = ma(0);
	m_p.translate(RPW/2.0+RPg+RPw/2.0);

	// the resonator configuration
	Circuit res = resdef();
	// vertical computational window
	Interval cwx(-CWx, CWx);
	// horizontal computational window
	Interval cwz(RPL/2.0-CWz, RPL/2.0+CWz);
	// QUEP setup + field discretization 
	QuepField qfld(res, Pol, cwx, NumMx, cwz, NumMz);

	// input: the fundamental mode in the lower left channel 
	qfld.input(LEFT, CC1, m_m);
	// calculate the light propagation through the resonator
	qfld.quepsim();

	fprintf(stderr, "\nScattered guided mode amplitudes:\n"); 
	Complex a;
	a = qfld.Aout(LEFT, m_m);
	fprintf(stderr, "R_m = %g + i %g,  |R_m|^2 = %g\n",
		a.re, a.im , a.sqabs());
	a = qfld.Aout(LEFT, m_p);
	fprintf(stderr, "R_p = %g + i %g,  |R_p|^2 = %g\n",
		a.re, a.im , a.sqabs());
	a = qfld.Aout(RIGHT, m_m);
	fprintf(stderr, "T_m = %g + i %g,  |T_m|^2 = %g\n",
		a.re, a.im , a.sqabs());
	a = qfld.Aout(RIGHT, m_p);
	fprintf(stderr, "T_p = %g + i %g,  |T_p|^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);
	fprintf(stderr, "Total scattered guided power:\n");
	fprintf(stderr, "P_back      = %g\n", pb);
	fprintf(stderr, "P_forw      = %g\n", pf);
	fprintf(stderr, "Sum P_out,g = %g\n\n", pb+pf);

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

	qfld.adjustphase();
	qfld.viewer(-DWx, DWx, DWz0, DWz1, 150, 150, 't', pc); 
	qfld.plot(fc, REP, -DWx, DWx, DWz0, DWz1, 150, 150, 't', pc); 
	qfld.plot(fc, IMP, -DWx, DWx, DWz0, DWz1, 150, 150, 't', pc); 
	qfld.plot(fc, MOD, -DWx, DWx, DWz0, DWz1, 150, 150, 't', pc);
	qfld.phasemap(fc, -DWx, DWx, DWz0, DWz1, 150, 150, 't', pc);
	qfld.movie(-DWx, DWx, DWz0, DWz1, 150, 150, 30, 't', pc);             
	qfld.fplot(fc, -DWx, DWx, DWz0, DWz1, 150, 150, 't', pc);
	qfld.fmovie(-DWx, DWx, DWz0, DWz1, 150, 150, 30, 't', pc);             

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