/* 
 * METRIC-Application:
 * A three layer waveguide, dispersion curves,
 * effective indices & mode angles versus film thickness
 */

#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   3.4    // film refractive index
#define Wgpnc   1.45   // cover refractive index
double  Wgpt;          // slab thickness 
#define Wavel   1.55   // vacuum wavelength 

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

#define Pmin    0.01   // parameter scan: minium, 
#define Pmax    7.0    //                 maximum level, 
#define NumP    699    //                 number of values

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

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 guided modes, depending on the varying film thickness, 
   store effective indices & mode angles vs. thickness in files */
int main()
{
	ModeArray ma;
	char name[20]    = "t___neff";
	char angname[20] = "t___ang";

	int omax = -1;
	name[1]    = polchr(Pol);
	angname[1] = polchr(Pol);

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

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

		// find the guided modes
		modeanalysis(wg, Pol, ma, 1);
		
		// store effective mode indices and mode angles
		for(int j=0; j<=ma.num-1; ++j)
		{
			int o = ma(j).nodes(); 
			if(o <= MaxOrd)
			{
				neff(o, pi) = ma(j).neff;	
				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); 
	for(int o=0; o<=omax; ++o)
	{
		int num = 0;
		name[2] = dig10(o);
		name[3] =  dig1(o);
		angname[2] = dig10(o);
		angname[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);
				++num;
			}
		}
		// ... effective indices 
		if(num >= 1)
		toxyf(name, oarg.subvector(num, 0), oval.subvector(num, 0));
		for(int pi=0; pi<=num-1; ++pi) 
			oval(pi) = acos(oval(pi)/Wgpnf)*180.0/PI;
		// ... and mode angles 
		if(num >= 1)
		toxyf(angname, oarg.subvector(num, 0), oval.subvector(num, 0));
	}
	return 0;
}