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

#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

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

/* 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, wavelength spectrum */
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;

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

	// 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;
	// remember the indices in the array of basis elements
	int top_f_eli, bot_b_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);
	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);
	apptoxyf("T", v, p);
	p = (fld(bot_b_eli).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");
        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);
	apptoxyf("pT", v, p);
tr.append(p);
        p = (fld.overlap_x(mb, BACK, CWz0)).sqabs(); tp += p;
        fprintf(stderr, " D: %g\n", p);
	apptoxyf("pD", v, p);
        fprintf(stderr, "Sum: %g\n", tp);
	apptoxyf("sumpP", v, tp);

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