/* 
 * METRIC-Application:
 * Wave propagation along a straight photonic crystal waveguide,
 * input & output via pieces of conventional waveguides,
 * rectangular lattice of square silicon rods in air
 */

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

#define Pol   TE      // light polarization
#define BPnr  3.4     // rod refractive index
#define BPnb  1.0     // background refractive index
#define BPr   0.150   // rod width /mum
#define BPp   0.600   // period of the photonic crystal /mum
#define BPw   0.500   // port waveguide width /mum
#define BPng  1.8     // port waveguide: core refractive index 
#define BPNr  4       // number of rods in the reflector regions around the 
                      // photonic crystal waveguide
#define BPNp  9       // number of rods along the pc waveguide (length)
#define Wavel 1.461   // vacuum wavelength 

#define NumMx 120     // vertical spectral discretization, number of terms 
#define NumMz 120     // horizontal spectral discretization, number of terms 
#define CW    1.2     // computational window: width of the border regions

#define DW    1.0     // display window
                      // (borders around the outer permittivity interfaces)
 
/* ------------------------------------------------------------------------ */

/* the photonic crystal configuration */
Circuit stdef()
{
	Circuit pcwg(BPNp-1+BPNp-1+1, 
	             BPNr+BPNr+BPNr-1+BPNr-1+2+2+1);
	int i=0; 
	double x = 0.0; pcwg.hx(i++) = x; 
	for(int j=1; j<=BPNp-1; ++j)
	{
		x += BPr; pcwg.hx(i++) = x; 
		x += BPp-BPr; pcwg.hx(i++) = x; 
	}
	x += BPr; pcwg.hx(i++) = x; 
	i=0; 
	double z = 0.0; pcwg.hz(i++) = z; 
	for(int j=1; j<=BPNr-1; ++j)
	{
		z += BPr; pcwg.hz(i++) = z; 
		z += BPp-BPr; pcwg.hz(i++) = z; 
	}
	z += BPr; pcwg.hz(i++) = z; 
	z += BPp-BPr/2.0-BPw/2.0; pcwg.hz(i++) = z; 
	z += BPw/2.0-BPr/2.0; pcwg.hz(i++) = z; 
	z += BPr; pcwg.hz(i++) = z; 
	z += BPw/2.0-BPr/2.0; pcwg.hz(i++) = z; 
	z += BPp-BPr/2.0-BPw/2.0; pcwg.hz(i++) = z; 
	z += BPr; pcwg.hz(i++) = z; 
	for(int j=1; j<=BPNr-1; ++j)
	{
		z += BPp-BPr; pcwg.hz(i++) = z; 
		z += BPr; pcwg.hz(i++) = z; 
	}
	pcwg.n.init(BPnb);
	int s, l, s0;
	s = (BPNr-1)*2+2+1;
	pcwg.n(0, s  ) = BPng; 
	pcwg.n(0, s+1) = BPng; 
	pcwg.n(0, s+2) = BPng; 
	pcwg.n(pcwg.nx+1, s  ) = BPng; 
	pcwg.n(pcwg.nx+1, s+1) = BPng; 
	pcwg.n(pcwg.nx+1, s+2) = BPng; 
	for(s=1; s<=2*BPNr-1; s+=2)
	{
		for(l=1; l<=2*BPNp-1; l+=2)
		{
			pcwg.n(l, s) = BPnr;
		}
	}
	s0 = 2*BPNr+4;
	for(s=s0+1; s<=s0+2*BPNr-1; s+=2)
	{
		for(l=1; l<=2*BPNp-1; l+=2)
		{
			pcwg.n(l, s) = BPnr;
		}
	}
	pcwg.lambda = Wavel;
	return pcwg;
}

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

/* QUEP simulation of the photonic crystal waveguide */ 
int main()
{
	
	fprintf(stderr, "\nPhotonic crystal waveguide:\n");	
	fprintf(stderr, "---------------------------\n");	
	fprintf(stderr, "lambda = %g mum\n", Wavel);
	fprintf(stderr, "n_b    = %g\n", BPnb);
	fprintf(stderr, "n_r    = %g\n", BPnr);
	fprintf(stderr, "p      = %g mum\n", BPp);
	fprintf(stderr, "r      = %g mum\n", BPr);
	fprintf(stderr, "n_g    = %g\n", BPng);
	fprintf(stderr, "w      = %g mum\n", BPw);
	fprintf(stderr, "Nr     = %d\n", BPNr);
	fprintf(stderr, "Np     = %d\n\n", BPNp);

	// the photonic crystal structure
	Circuit pcwg = stdef();
	// a plot of the refractive index profile
	pcwg.plot('0', '0');

	// computational window
        Interval cw(pcwg.hz(0)-CW, pcwg.hz(pcwg.nz)+CW);
	// QUEP setup + spectral discretization
	QuepField qfld(pcwg, Pol, cw, NumMx, cw, NumMz);

        // input: the fundamental mode in the top channel 
	qfld.input(TOP, CC1, 0);
        // calculate the light propagation through the configuration
        qfld.quepsim();

	Complex a;
	fprintf(stderr, "\nScattered fundamental mode amplitudes:\n");
	a = qfld.Aout(BOTTOM, 0);
	fprintf(stderr, "T = %g + i %g,  |T|^2 = %g\n", a.re, a.im, a.sqabs());
	a = qfld.Aout(TOP, 0);
	fprintf(stderr, "R = %g + i %g,  |R|^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);

        pd = qfld.Pgout(BOTTOM);
        pu = qfld.Pgout(TOP);
        fprintf(stderr, "Total scattered guided power:\n");
        fprintf(stderr, "P_down      = %g\n", pd);
        fprintf(stderr, "P_up        = %g\n", pu);
        fprintf(stderr, "Sum P_out,g = %g\n\n", pd+pu);

	double dwl, dwr, dwb, dwt;
	dwl = pcwg.hz(0)-DW; 
	dwr = pcwg.hz(pcwg.nz)+DW;
	dwb = dwl; 
	dwt = dwr;
	Fcomp fc = principalcomp(Pol);
	char pc = polchr(Pol);

	qfld.adjustphase();
        qfld.viewer(dwb, dwt, dwl, dwr, 250, 250, 't', pc);
        qfld.plot(fc, REP, dwb, dwt, dwl, dwr, 250, 250, 't', pc);
        qfld.plot(fc, IMP, dwb, dwt, dwl, dwr, 250, 250, 't', pc);
        qfld.plot(fc, MOD, dwb, dwt, dwl, dwr, 250, 250, 't', pc);
        qfld.phasemap(fc, dwb, dwt, dwl, dwr, 250, 250, 't', pc);
        qfld.movie(dwb, dwt, dwl, dwr, 250, 250, 30, 't', pc);
        qfld.fplot(fc, dwb, dwt, dwl, dwr, 250, 250, 't', pc);
        qfld.fmovie(dwb, dwt, dwl, dwr, 250, 250, 30, 't', pc);

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