/* 
 * METRIC-Application:
 * Crossing between perpendicular waveguides
 */

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

#define Pol    TE     // light polarization 
#define CPnb   1.45   // background refractive index 
#define CPng   3.40   // waveguide cores: refractive index 
#define CPw    0.20   // width of the horizontal waveguide /mum 
#define CPv    0.20   // width of the vertical waveguide /mum 
#define Wavel  1.55   // vacuum wavelength /mum

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

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

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

/* a three layer slab waveguide of thickness t */
Waveguide wgdef(double t)
{
	Waveguide g(1);
	g.hx(0) = -t/2.0;
	g.hx(1) =  t/2.0;
	g.n(0) = CPnb;
	g.n(1) = CPng;
	g.n(2) = CPnb;
	g.lambda = Wavel;
	return g;
}

/* the waveguide crossing */
Circuit circuitdef()
{
	Circuit cr;
	cr = Circuit(1, 1);
	cr.hx(0) = -CPw/2.0;
	cr.hx(1) =  CPw/2.0;
	cr.hz(0) = -CPv/2.0;
	cr.hz(1) =  CPv/2.0;
	cr.n(0, 0) = CPnb; cr.n(0, 1) = CPng; cr.n(0, 2) = CPnb;
	cr.n(1, 0) = CPng; cr.n(1, 1) = CPng; cr.n(1, 2) = CPng;
	cr.n(2, 0) = CPnb; cr.n(2, 1) = CPng; cr.n(2, 2) = CPnb;
	cr.lambda = Wavel;
	return cr;
}

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

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

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

/* crossing simulation */
int main()
{
	fprintf(stderr, "\nWaveguide crossing:\n");	
	fprintf(stderr, "-------------------\n");	
	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", CPnb);
	fprintf(stderr, "n_g    = %g\n", CPng);
	fprintf(stderr, "w      = %g mum\n", CPw);
	fprintf(stderr, "v      = %g mum\n", CPv);

	fprintf(stderr, "\nHorizontal waveguide:\n");
	ModeArray mh;
	Waveguide hwg = wgdef(CPw);
	modeanalysis(hwg, Pol, mh);
	if(mh.num <= 0) crosserror("Horizontal waveguide: No mode found");
	fprintf(stderr, "Vertical waveguide:\n");
	ModeArray mv;
	Waveguide vwg = wgdef(CPv);
	modeanalysis(vwg, Pol, mv);
	if(mv.num <= 0) crosserror("Vertical waveguide: No mode found");
	fprintf(stderr, "\n");	

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

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

	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n"); 
	Complex a; 
	double p;
	a = qfld.Aout(LEFT, 0);
	p = qfld.Pout(LEFT, 0);
	fprintf(stderr, "B = %g + i %g,  |B|^2 = %g\n", a.re, a.im, p);
	a = qfld.Aout(RIGHT, 0);
	p = qfld.Pout(RIGHT, 0);
	fprintf(stderr, "F = %g + i %g,  |F|^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);
	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);
	
	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);
	pd = qfld.Pgout(BOTTOM);
	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_down      = %g\n", pd);
	fprintf(stderr, "P_up        = %g\n", pu);
	fprintf(stderr, "Sum P_out,g = %g\n\n", pb+pf+pd+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; 
}