/* 
 * METRIC-Application:
 * Three coupled circular microcavities, triangular arrangement, 
 * HCMT model, WGM basis, supermode computation
 */

#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.5621;   // target 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

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

/* 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, supermode analysis */
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(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);

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

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

	// HCMT scheme, supermode analysis
        Cvector om;
        om = fld.supermodes();

        fprintf(stderr, "\n");
        for(int j=0; j<=om.nel-1; ++j)
        {
		Complex o = fld.supermode(j);
		double l = val_lambda(o.re);
		fprintf(stderr, "omega_%d = %g + i %g,  lambda_r = %g,  Q = %g,  FWHM_lambda = %g\n", j, o.re, o.im, l, theQ(o), FWHM_lambda(o));
                Fcomp pc = principalcomp(Pol);
		fld.adjustphase(CWx0, CWx1, CWz0, CWz1, 350, 350);
                fld.plot(pc, MOD, CWx0, CWx1, CWz0, CWz1, 350, 350, dig10(j), dig1(j));
                fld.plot(pc, REP, CWx0, CWx1, CWz0, CWz1, 350, 350, dig10(j), dig1(j));
                fld.viewer(CWx0, CWx1, CWz0, CWz1, 350, 350, dig10(j), dig1(j));
        }
        fprintf(stderr, "\nOk.\n");
        return 0;
}