/* 
 * METRIC-Application:
 * A circular microring-resonator, two parallel bus waveguides, symmetric setting,
 * modeling by hybrid analytical / numerical coupled mode theory
 * wavelength spectrum, fast evaluation by interpolation of HCMT operators
 */

#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 
double Wavel = 1.55;   // vacuum wavelength /mum

// 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)


#define Pmin   1.520               // parameter scan: minimum,
#define Pmax   1.620               //                 maximum,
#define H_0    (Pmax-Pmin)/200.0   // initial scan stepsize
#define H_min  (Pmax-Pmin)/1000.0  // minimum local scan stepsize

               // interpolated parameter scan:
#define Nv0    (Pmin+(Pmax-Pmin)*0.25) // nodal value 0  
#define Nv1    (Pmin+(Pmax-Pmin)*0.75) // nodal value 1  

// remember the indices in the array of basis elements 
int tEli;
int bEli;

// cavity bend mode, integer angular wavenumbers, uniform for all wavelengths 
int KappaR = 0;

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

/* 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;
}

HcmtField setup()
{
	// 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, "mrspec: 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, "mrspec: 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;

	// 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; 
	tEli = 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; 
	bEli = fld.addcf(el);

	// cavity, clockwise propagating bend mode 
	if(KappaR == 0)
                KappaR = (int)(floor(bma(0).beta*bma(0).R+0.5));

	el = HcmtElement(CCC, bma(0), 0.0, 0.0, FEMstep/PR, KappaR);
	fld.addcf(el);

	return fld;
}

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

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

	HcmtField fld0, fld1, fld;
        Cmatrix M0, M1, M; 

	// field containers for interpolation procedure
        Wavel = 0.5*(Nv0+Nv1);
        fld = setup();
        fld.init();
	// ... at nodal value 0
        Wavel = Nv0;
        fld0 = setup();
        fld0.init();
        fld0.setupsystem(M0);
	// ... at nodal value 1
        Wavel = Nv1;
        fld1 = setup();
        fld1.init();
        fld1.setupsystem(M1);

// wavelength scan: 
// rough initial survey, locating transmittance minima through bisection 
Dvector wl, tr;
wl.strip();
tr.strip();
int done=0;
while(done==0)
{
done = next(Interval(Pmin, Pmax), H_0, H_min, wl, tr);
if(done==0)
{
	double v = wl(wl.nel-1);
        Wavel = v;
        fprintf(stderr, "\n++++  lambda = %g\n", Wavel);

	// fast spectrum evaluation: interpolate HCMT operators
	M = add(M0, mult(sub(M1, M0), (Wavel-Nv0)/(Nv1-Nv0)));
        fld.solve(M);

	// result: output power
	fprintf(stderr, "\n");
	double tp = 0.0, p;
	fprintf(stderr, "Squared output amplitudes: \n");
	p = (fld(tEli).aout).sqabs(); tp += p;
	fprintf(stderr, " T: %g\n", p);
	apptoxyf("T", v, p);
tr.append(p);
	p = (fld(bEli).aout).sqabs(); tp += p;
	fprintf(stderr, " D: %g\n", p);
	apptoxyf("D", v, p);
	fprintf(stderr, "Sum: %g\n", tp);
	apptoxyf("sumP", v, tp);

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