/* 
 * METRIC-Application:
 * A circular microring-resonator, two parallel bus waveguides, symmetric setting,
 * modeling by hybrid analytical / numerical coupled mode theory
 */

#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 /mum 
#define PR    7.5    // cavity radius, outer rim /mum 
#define Pd    0.75   // cavity core, width /mum 
#define Wavel 1.5621 // vacuum wavelength /mum

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


// HCMT basis fields, parameters of the FEM discretization of amplitudes
#define FEMstep 0.1
#define Z0 -10.0
#define ZN  10.0

// evaluation of HCMT integrals: computational window
#define CWz0   -10.0   
#define CWz1    10.0 
#define CWx0   -10.0 
#define CWx1    10.0 

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

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

/* 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 cavity core */
Waveguide bwgdef()
{
	Waveguide b(1);
	b.hx(1) =  0.0;  
	b.hx(0) =  -Pd;  
	b.n(2) = Pnb;
	b.n(1) = Png;
	b.n(0) = Pnb;
	b.lambda = Wavel;
	return b;
}

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

/* microresonator simulation */
int main()
{
	// crop mode profiles where rel. level is below MCROP
	MCROP = 1.0e-6;
	// disregard bend mode profiles outside radial interval
	BM_R_CROP_IN  = PR-3.0;
	BM_R_CROP_OUT = PR+10.0;

	// the composite resonator structure
	OvlStruct str(Pnb*Pnb, Wavel);
	Overlay o;
	o = Overlay(HLAY, Png, Pw,  PR+Pg+0.5*Pw);
	str.place(o);
	o = Overlay(HLAY, Png, Pw, -PR-Pg-0.5*Pw);
	str.place(o);
	o = Overlay(RING, Png, PR, Pd, 0.0, 0.0);
	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 ma;
	modeanalysis(wg, Pol, ma);
	if(ma.num != 1) { fprintf(stderr, "mrres: bus cores are not singlemode.\n"); exit(1); }

	// bend mode of the cavity
	Waveguide bwg = bwgdef();
	BDModeArray bma;
	bdmsolve(bwg, PR, 0.0, 0.0, Pol, 1, 20, 10, bma, 1);
	if(bma.num != 1) { fprintf(stderr, "mrres: cavity core is not singlemode.\n"); exit(1); }
	bma(0).discretize(BM_R_CROP_IN, BM_R_CROP_OUT, 5000);

	// add basis fields: 
	HcmtElement el;
	// remember the indices in the array of basis elements
	int top_f_eli, bot_b_eli, cav_c_eli; 

	// upper bus waveguide, forward mode, unit input amplitude 
	Mode mt = ma(0); 
	mt.translate( PR+Pg+Pw/2.0);
	el = HcmtElement(MHF, mt, Z0, ZN, FEMstep);
	el.ain = CC1; 
	top_f_eli = fld.addcf(el);

	// lower bus waveguide, backward mode, zero input 
	Mode mb = ma(0); 
	mb.translate(-PR-Pg-Pw/2.0);
	el = HcmtElement(MHB, mb, Z0, ZN, FEMstep);
	el.ain = CC0; 
	bot_b_eli = fld.addcf(el);

	// cavity, clockwise propagating bend mode 
	el = HcmtElement(CCC, bma(0), 0.0, 0.0, FEMstep/PR);
	cav_c_eli = fld.addcf(el);

	// solve the problem according to the HCMT scheme
	fld.solve();

	// result: output power
	fprintf(stderr, "\n");
	double tp = 0.0, p;
	fprintf(stderr, "Squared output amplitudes: \n");
	p = (fld(top_f_eli).aout).sqabs(); tp += p;
	fprintf(stderr, " T: %g\n", p);
	p = (fld(bot_b_eli).aout).sqabs(); tp += p;
	fprintf(stderr, " D: %g\n", p);
	fprintf(stderr, "Sum: %g\n", tp);

	fprintf(stderr, "\n");
	tp = 0.0;
	fprintf(stderr, "Squared output amplitudes, projected: \n");
	p = (fld.overlap_x(mt, FORW, CWz1)).sqabs(); tp += p;
	fprintf(stderr, " T: %g\n", p);
	p = (fld.overlap_x(mb, BACK, CWz0)).sqabs(); tp += p;
	fprintf(stderr, " D: %g\n", p);
	fprintf(stderr, "Sum: %g\n", tp);

	// amplitude functions -> data files
	fprintf(stderr, "\n");
	fld(top_f_eli).amp2files('t', 'f');
	fld(bot_b_eli).amp2files('b', 'b');
	fld(cav_c_eli).amp2files('c', 'c');

	// field plot 
	fprintf(stderr, "\n");
	Fcomp pc = principalcomp(Pol);
	fld.plot(pc, REP, -DWx, DWx, -DWz, DWz, 350, 350, '0', '0');
	fld.plot(pc, MOD, -DWx, DWx, -DWz, DWz, 350, 350, '0', '0');
	fld.viewer(-DWx, DWx, -DWz, DWz, 350, 350, '0', '0');

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