/* 
 * METRIC-Application:
 * Three coupled circular microcavities, triangular arrangement, 
 * one bus waveguide, symmetric setting
 * HCMT model, WGM basis
 */

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

#define Pol    TE
#define Pnb   1.0    // background refractive index
#define Png   1.5    // cores: refractive index
#define Pw    0.6    // bus waveguides, width /mum
#define Pg    0.3    // gaps wg-cav /mum
#define Ps    0.3    // gap cav-cav /mum
#define PR    7.5    // cavity radii, outer rim /mum
#define Pd    0.75   // cavity cores, width /mum
double Wavel = 1.55969;  // vacuum wavelength /mum

// #define MinWl 1.48   // include WGMs in this wavelength interval
// #define MaxWl 1.64

// evaluation of HCMT integrals: computational window
#define Mar     (10.0-PR)
#define CWz0   -Mar
#define CWz1    (2.0*Pg+4.0*PR+Mar)
#define CWx0   -(2.0*PR+Ps/2.0+Mar)
#define CWx1    (2.0*PR+Ps/2.0+Mar)

#define Nispwl 4   // number of integration steps (10 pt Gaussian quadrature)
                   // per wavelength

// HCMT basis fields, parameters of the FEM discretization of amplitudes
#define FEMstep 0.2
#define X0  CWx0
#define XN  CWx1

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

/* prototype for the bus waveguides */
Waveguide wgdef()
{
        Waveguide g(1);
        g.hx(0) = -Pw/2.0;
        g.hx(1) =  Pw/2.0;
        g.n(0) = Pnb;
        g.n(1) = Png;
        g.n(2) = Pnb;
        g.lambda = Wavel;
        return g;
}

/* the core of the cavity rings */
Waveguide rwgdef()
{
        Waveguide rwg(1);
        rwg.hx(0) =  -Pd;
        rwg.hx(1) =  0.0;
        rwg.n(0) = Pnb;
        rwg.n(1) = Png;
        rwg.n(2) = Pnb;
        rwg.lambda = Wavel;

        return rwg;
}

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

/* microresonator simulation, scattering problem */
int main()
{
	MCROP = 1.0e-6;
	WM_R_CROP_IN  = 0.0;
	WM_R_CROP_OUT = PR*4.0;

	// compute the cavity positions 
	double alpha;
	alpha = asin((PR+Ps/2.0)/(2.0*PR+Pg));
	double c0x0, c0z0, c1x0, c1z0, c2x0, c2z0;
	c0x0 = 0.0;
	c0z0 = Pg+PR;
	c1x0 =   Ps/2.0+PR;
	c1z0 = Pg+PR+(2.0*PR+Pg)*cos(alpha);
	c2x0 = -(Ps/2.0+PR);
	c2z0 = Pg+PR+(2.0*PR+Pg)*cos(alpha);

	// the composite circuit
	OvlStruct str(Pnb*Pnb, Wavel);
	Overlay o;
	o = Overlay(VLAY, Png, Pw,  -0.5*Pw);
	str.place(o);
	o = Overlay(RING, Png, PR, Pd, c0x0, c0z0);
	str.place(o);
	o = Overlay(RING, Png, PR, Pd, c1x0, c1z0);
	str.place(o);
	o = Overlay(RING, Png, PR, Pd, c2x0, c2z0);
	str.place(o);

	// HCMT field container, initialization
	HcmtField fld(str, Pol, Wavel, 
	              Interval(CWx0, CWx1), Interval(CWz0, CWz1), Nispwl);

	// mode of the bus core
	Waveguide wg = wgdef();
	ModeArray mh;
	modeanalysis(wg, Pol, mh);
	mh(0).translate(-Pw/2.0);

        // WGM(s) of the cavity 
        WGModeArray wgm;
        Waveguide rwg = rwgdef();
        wgmsolve(rwg, PR, 0.0, 0.0, Pol, Wavel, 1, -1, -1, wgm, 1);
	while(wgm.num > 1) wgm.remove(wgm.num-1); // keep only one WGM
//      wgmsolve(rwg, PR, 0.0, 0.0, Pol, Interval(MinWl, MaxWl), 1, -1, -1, wgm, 1);
	WGModeArray wgmr;
	wgmr = wgm;
	for(int j=0; j<=wgmr.num-1; ++j) 
	{
		wgmr(j).reverse();
		wgm.add(wgmr(j));
	}
	for(int j=0; j<=wgm.num-1; ++j) wgm(j).discretize(0.0, 4.0*PR, 4.0*PR/(Pd/20.0));

	// add HCMT basis fields
	HcmtElement el;

	// upward travelling bus mode
	el = HcmtElement(MVU, mh(0), X0, XN, FEMstep);
	el.ain = CC1; 
	int wg_u_eli = fld.addcf(el);

	// downward travelling bus mode
	el = HcmtElement(MVD, mh(0), X0, XN, FEMstep);
	el.ain = CC0; 
	int wg_d_eli = fld.addcf(el);

	// WGMs of the cavities
	wgmr = wgm;
	for(int j=0; j<=wgm.num-1; ++j)
	{ 
		wgmr(j).translate(c0x0, c0z0);
		el = HcmtElement(WGM, wgmr(j));
		fld.addcf(el);
	}
	wgmr = wgm;
	for(int j=0; j<=wgm.num-1; ++j)
	{ 
		wgmr(j).translate(c1x0, c1z0);
		el = HcmtElement(WGM, wgmr(j));
		fld.addcf(el);
	}
	wgmr = wgm;
	for(int j=0; j<=wgm.num-1; ++j)
	{ 
		wgmr(j).translate(c2x0, c2z0);
		el = HcmtElement(WGM, wgmr(j));
		fld.addcf(el);
	}

	// solve the scattering problem
	fld.solve();

	fprintf(stderr, "\n");
	double tp = 0.0, p;
	fprintf(stderr, "Squared output amplitudes: \n");
	p = (fld(wg_u_eli).aout).sqabs(); tp += p;
	fprintf(stderr, " T: %g\n", p);
	p = (fld(wg_d_eli).aout).sqabs(); tp += p;
	fprintf(stderr, " R: %g\n", p);
	fprintf(stderr, "Sum: %g\n", tp);

	Fcomp pc = principalcomp(Pol);
	fld.adjustphase(CWx0, CWx1, CWz0, CWz1, 350, 350);
	fld.plot(pc, MOD, CWx0, CWx1, CWz0, CWz1, 350, 350, '0', '0');
	fld.plot(pc, REP, CWx0, CWx1, CWz0, CWz1, 350, 350, '0', '0');
        fld.viewer(CWx0, CWx1, CWz0, CWz1, 350, 350, '0', '0');

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