/* 
 * METRIC-Application:
 * A three layer waveguide, wavelength dependence of propagation constants 
 * and effective mode indices, rigorous curves & perturbational evaluation 
 * of wavelength derivatives 
 */

#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
#define Wgpt    0.5     // slab thickness /mum
double  Wavel = 1.55;   // vacuum wavelength /mum

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

#define Pbeg    0.4     // Parameter variations: minimum value, 
#define Pend    1.6     //                       maximum value, 
#define NumP    240     //                       number of wavelengths

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

/* 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 propagation constants and effective mode indices 
   versus the vacuum wavelength parameter */
int main()
{
	ModeArray ma;
	char propname[20] = "t___beta";
	char neffname[20] = "t___neff";
	char pcgrname[20] = "t_dbdl";
	char nfgrname[20] = "t_dndl";

	propname[1] = polchr(Pol);
	neffname[1] = polchr(Pol);
	pcgrname[1] = polchr(Pol);
	nfgrname[1] = polchr(Pol);

	fprintf(stderr, 
	"\nWavelength gradients of propagation constants  >> %s.xyf\n", 
	pcgrname);
	fprintf(stderr, 
	"Wavelength gradients of effective mode indices >> %s.xyf\n", 
	nfgrname);

	Dmatrix beta(MaxOrd+1, NumP+1); 
	Dmatrix neff(MaxOrd+1, NumP+1); 
	Dvector pval(NumP+1); 
	int omax = -1;
	beta.init(0.0);
	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);
		Wavel = pval(pi);	

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

		// find the guided modes 
		modeanalysis(wg, Pol, ma, 1);
		
		// store propagation constants and effective mode indices
		for(int j=0; j<=ma.num-1; ++j)
		{
			int o = ma(j).nodes(); 
			if(o <= MaxOrd)
			{
				beta(o, pi) = ma(j).beta;	
				neff(o, pi) = ma(j).neff;	
				if(o > omax) omax = o;
				// at about 10 intermediate wavelengths ... 
				if(pi%(NumP/10)-(NumP/10)/2 == 0)
				{  
					// ... evaluate the wavel. derivatives
					//     of the propagation constant 
					//     and of the effective mode index
					Complex db, dn;
					db = ma(j).phaseshift(1.0, FORW);
					dn = ma(j).neffshift(1.0, FORW);
					// ... and append them to files
					apptoxyf(pcgrname, 
					         Wavel, ma(j).beta, db.re);
					apptoxyf(nfgrname, 
					         Wavel, ma(j).neff, dn.re);
				}
			}
		}
		int s = NumP/60;
                if(pi % s == 0) fprintf(stderr, ".");
	}
	fprintf(stderr, "\n");
	// save the rigorous data to files, individually per mode order, for ...
	Dvector oarg(NumP+1); 
	Dvector bval(NumP+1); 
	Dvector nval(NumP+1); 
	for(int o=0; o<=omax; ++o)
	{
		int num = 0;
		propname[2] = dig10(o);
		propname[3] =  dig1(o);
		neffname[2] = dig10(o);
		neffname[3] =  dig1(o);
		for(int pi=0; pi<=NumP; ++pi)
		{
			if(neff(o, pi) > 0.0)	
			{
				oarg(num) = pval(pi);
				bval(num) = beta(o, pi);
				nval(num) = neff(o, pi);
				++num;
			}
		}
		if(num >= 1)
		{
		// ... propagation constants
		toxyf(propname, oarg.subvector(num, 0), bval.subvector(num, 0));
		// ... and effective mode indices	
		toxyf(neffname, oarg.subvector(num, 0), nval.subvector(num, 0));
		}
	}
	fprintf(stderr, "Ok.\n");
	return 0;
}