/* 
 * METRIC-Application:
 * Wave propagation along a 90^o bend of a photonic crystal waveguide,
 * input and 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 BPN   4       // number of rods in the reflector regions around the 
                      // photonic crystal waveguide
#define Wavel 1.461   // vacuum wavelength 

                      // the discretized radiation field
#define NumM  120     // number of expansion terms for both directions
#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 pcb(BPN+BPN+BPN-1+BPN-1+2+2+1, 
	            BPN+BPN+BPN-1+BPN-1+2+2+1);
	int i=0; 
	double x = 0.0; pcb.hx(i++) = x; 
	for(int j=1; j<=BPN-1; ++j)
	{
		x += BPr; pcb.hx(i++) = x; 
		x += BPp-BPr; pcb.hx(i++) = x; 
	}
	x += BPr; pcb.hx(i++) = x; 
	x += BPp-BPr/2.0-BPw/2.0; pcb.hx(i++) = x; 
	x += BPw/2.0-BPr/2.0; pcb.hx(i++) = x; 
	x += BPr; pcb.hx(i++) = x; 
	x += BPw/2.0-BPr/2.0; pcb.hx(i++) = x; 
	x += BPp-BPr/2.0-BPw/2.0; pcb.hx(i++) = x; 
	x += BPr; pcb.hx(i++) = x; 
	for(int j=1; j<=BPN-1; ++j)
	{
		x += BPp-BPr; pcb.hx(i++) = x; 
		x += BPr; pcb.hx(i++) = x; 
	}
	pcb.hz = pcb.hx;
	pcb.n.init(BPnb);
	int s, l, l0, s0;
	s = (BPN-1)*2+2+1;
	pcb.n(0, s  ) = BPng; 
	pcb.n(0, s+1) = BPng; 
	pcb.n(0, s+2) = BPng; 
	l = (BPN-1)*2+2+1;
	pcb.n(l,   pcb.nz+1) = BPng; 
	pcb.n(l+1, pcb.nz+1) = BPng; 
	pcb.n(l+2, pcb.nz+1) = BPng; 
	for(s=1; s<=2*BPN-1; s+=2)
	{
		for(l=1; l<=2*BPN-1; l+=2)
		{
			pcb.n(l, s) = BPnr;
		}
		l = 2*BPN+2; pcb.n(l, s) = BPnr;
		l0 = 2*BPN+4;
		for(l=l0+1; l<=l0+2*BPN-1; l+=2)
		{
			pcb.n(l, s) = BPnr;
		}
	}
	s = BPN*2+2;
	l0 = 2*BPN+4;
	for(l=l0+1; l<=l0+2*BPN-1; l+=2)
	{
		pcb.n(l, s) = BPnr;
	}
	s0 = 2*BPN+4;
	for(s=s0+1; s<=s0+2*BPN-1; s+=2)
	{
		for(l=1; l<=2*BPN-1; l+=2)
		{
			pcb.n(l, s) = BPnr;
		}
		l0 = 2*BPN+4;
		for(l=l0+1; l<=l0+2*BPN-1; l+=2)
		{
			pcb.n(l, s) = BPnr;
		}
	}
	pcb.lambda = Wavel;
	return pcb;
}

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

/* simulation of the photonic crystal bend */ 
int main()
{
	
	fprintf(stderr, "\nPhotonic crystal bend:\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, "N      = %d\n\n", BPN);

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

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

        // input: the fundamental mode in the bottom channel
	qfld.input(BOTTOM, CC1, 0);
        // calculate the light propagation along the bend 
        qfld.quepsim();

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

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

	double dwl, dwr, dwb, dwt;
	dwl = pcb.hz(0)-DW; 
	dwr = pcb.hz(pcb.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;
}