/* 
 * METRIC-Application:
 * A three layer waveguide, attenuating core,
 * dispersion curves, perturbational computation of mode attenuation 
 */

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

#define Pol     TE      // light polarization
#define Wgpns   1.45    // substrate refractive index
#define Wgpnf   1.98    // film refractive index
#define Wgpnc   1.0     // cover: air
double  Wgpt  = 0.5;    // slab thickness /mum
#define Wavel   1.55    // vacuum wavelength /mum
#define WgpAtt  1.0e-4  // attenuation constant of the core, 
                        // intensity attenuation (in 1/mum = 10^4/cm) 
                        // of a plane wave in the bulk material

#define MaxOrd 100      // highest order of a recorded mode 

#define Pbeg 0.01       // Parameter variations: minimum value, 
#define Pend 2.0        //                       maximum value, 
#define NumP 200        //                       number of points

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

/* waveguide definition */
Waveguide wgdef()
{
	Waveguide g(1);
	g.hx(0) = 0.0;
	g.hx(1) = Wgpt;
	g.n(0) = Wgpns;
	g.n(1) = Wgpnf;
	g.n(2) = Wgpnc;
	g.lambda = Wavel;
	return g;
}

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

/* calculate effective indices and relative attenuation levels
   (intensity attenuation of the guided modes relative to the intensity 
    attenuation of plane waves in the bulk core material) 
   versus the layer thickness parameter */
int main()
{
	ModeArray ma;
	char name[20] = "t___neff";
	char attname[20] = "t___att";

	name[1] = polchr(Pol);
	attname[1] = polchr(Pol);

	Dmatrix neff(MaxOrd+1, NumP+1); 
	Dmatrix att(MaxOrd+1, NumP+1); 
	Dvector pval(NumP+1); 
	int omax = -1;
	neff.init(0.0);
	
	fprintf(stderr, "\nDispersion curves:\n");
	for(int pi=0; pi<=NumP; ++pi)
	{
		pval(pi) = Pbeg+((double)pi)/((double)NumP)*(Pend-Pbeg);
		Wgpt = pval(pi);	

		// define the waveguide
		Waveguide wg = wgdef();

		// the permittivity perturbation, loss in the core layer
		Perturbation pert = attenuation(WgpAtt, wg.n(1), wg.lambda, 
		                                wg.layer(1));

		// find the guided modes 
		modeanalysis(wg, Pol, ma, 1);
		
		// store effective mode indices, evaluate the perturbation
		for(int j=0; j<=ma.num-1; ++j)
		{
			int o = ma(j).nodes(); 
			if(o <= MaxOrd)
			{
				neff(o, pi) = ma(j).neff;	
				// mode intensity attenuation,
				// relative to the bulk attenuation 
				att(o, pi) = -ma(j).phaseshift(pert, FORW).im
				              *2.0/WgpAtt; 
				if(o > omax) omax = o;
			}
		}
		int s = NumP/60;
                if(pi % s == 0) fprintf(stderr, ".");
	}
	fprintf(stderr, "\n");
	// write the data to files, individually per mode order, for ...
	Dvector oarg(NumP+1); 
	Dvector oval(NumP+1); 
	Dvector aval(NumP+1); 
	for(int o=0; o<=omax; ++o)
	{
		int num = 0;
		name[2] = dig10(o);
		name[3] =  dig1(o);
		attname[2] = dig10(o);
		attname[3] =  dig1(o);
		for(int pi=0; pi<=NumP; ++pi)
		{
			if(neff(o, pi) > 0.0)	
			{
				oarg(num) = pval(pi);
				oval(num) = neff(o, pi);
				aval(num) = att(o, pi);
				++num;
			}
		}
		if(num >= 1)
		{
		// ... effective mode indices, 
		toxyf(name, oarg.subvector(num, 0), oval.subvector(num, 0));
		// ... and attenuation levels
		toxyf(attname, oarg.subvector(num, 0), aval.subvector(num, 0));
		}
	}
	fprintf(stderr, "Ok.\n");
	return 0;
}