/*
 * WMM - Wave-matching method for mode analysis of dielectric waveguides
 * https://wmm.computational-photonics.eu/
 */

/*
 * wmmgam.cpp
 * spectral discretization 
 */

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"waveg.h"
#include"maxweq.h"
#include"wmmtrif.h"
#include"wmmgam.h"

#define PI 3.14159265358979323846

double GPxa, GPxr, GPxb;
double GPya, GPyr, GPyb;
double GPw;

double GFDiff    = 0.01;
double MaxWvpA   = 3.0;

int Vlsg;

#define DeltaAlpha 0.01
#define DeltaF 0.01

#define AlphaStart PI/4.0;
#define FStart 1.0;

double gamalph(Fstype typ, double pvz, double qvz, double a, double da) 
{
	double p1, p2, q1, q2;
	double g, n1, n2;
	Rect r;

	p1 = pvz*GPw*cos(a-0.5*da);
	q1 = qvz*GPw*sin(a-0.5*da);
	p2 = pvz*GPw*cos(a+0.5*da);
	q2 = qvz*GPw*sin(a+0.5*da);

	Trifun tf1(0, 1.0, typ, 1.0, '-', p1, 1, q1, 1, GPxr, GPyr, '+', 1.0);
	Trifun tf2(0, 1.0, typ, 1.0, '-', p2, 1, q2, 1, GPxr, GPyr, '+', 1.0);
	r.x0 = GPxa;
	r.y0 = GPya;
	r.x1 = GPxb;
	r.y1 = GPyb;

	n1 = tfnrm(tf1, FLD, tf1, FLD, r);
	n1 = sqrt(n1);
	n2 = tfnrm(tf2, FLD, tf2, FLD, r);
	n2 = sqrt(n2);
	g = tfnrm(tf1, FLD, tf2, FLD, r)/n1/n2;
	return fabs(g);
}

double ddagamalph(Fstype typ, double pvz, double qvz, double a)
{
	double gp, g0, gm, ddag;
	gp = gamalph(typ, pvz, qvz, a, DeltaAlpha);
	g0 = gamalph(typ, pvz, qvz, a, 0);
	gm = gamalph(typ, pvz, qvz, a, -DeltaAlpha);
	ddag = (gp-2.0*g0+gm)/DeltaAlpha/DeltaAlpha;
	return ddag;
}


double gamfpf(Fstype typ, double pvz, double qvz, double f, double df) 
{
	double p1, p2, q1, q2;
	double g, n1, n2;
	Rect r;

	p1 = pvz*GPw*sinh(f-0.5*df);
	q1 = qvz*GPw*cosh(f-0.5*df);
	p2 = pvz*GPw*sinh(f+0.5*df);
	q2 = qvz*GPw*cosh(f+0.5*df);

	Trifun tf1(0, 1.0, typ, 1.0, '-', p1, 1, q1, 1, GPxr, GPyr, '+', 1.0);
	Trifun tf2(0, 1.0, typ, 1.0, '-', p2, 1, q2, 1, GPxr, GPyr, '+', 1.0);
	r.x0 = GPxa;
	r.y0 = GPya;
	r.x1 = GPxb;
	r.y1 = GPyb;

	n1 = tfnrm(tf1, FLD, tf1, FLD, r);
	n1 = sqrt(n1);
	n2 = tfnrm(tf2, FLD, tf2, FLD, r);
	n2 = sqrt(n2);
	g = tfnrm(tf1, FLD, tf2, FLD, r)/n1/n2;
	return fabs(g);
}

double ddfgamfpf(Fstype typ, double pvz, double qvz, double f)
{
	double gp, g0, gm, ddfg;
	gp = gamfpf(typ, pvz, qvz, f, DeltaF);
	g0 = gamfpf(typ, pvz, qvz, f, 0);
	gm = gamfpf(typ, pvz, qvz, f, -DeltaF);
	ddfg = (gp-2.0*g0+gm)/DeltaF/DeltaF;
	return ddfg;
}


double gamfqf(Fstype typ, double pvz, double qvz, double f, double df) 
{
	double p1, p2, q1, q2;
	double g, n1, n2;
	Rect r;

	p1 = pvz*GPw*cosh(f-0.5*df);
	q1 = qvz*GPw*sinh(f-0.5*df);
	p2 = pvz*GPw*cosh(f+0.5*df);
	q2 = qvz*GPw*sinh(f+0.5*df);

	Trifun tf1(0, 1.0, typ, 1.0, '-', p1, 1, q1, 1, GPxr, GPyr, '+', 1.0);
	Trifun tf2(0, 1.0, typ, 1.0, '-', p2, 1, q2, 1, GPxr, GPyr, '+', 1.0);
	r.x0 = GPxa;
	r.y0 = GPya;
	r.x1 = GPxb;
	r.y1 = GPyb;

	n1 = tfnrm(tf1, FLD, tf1, FLD, r);
	n1 = sqrt(n1);
	n2 = tfnrm(tf2, FLD, tf2, FLD, r);
	n2 = sqrt(n2);
	g = tfnrm(tf1, FLD, tf2, FLD, r)/n1/n2;
	return fabs(g);
}

double ddfgamfqf(Fstype typ, double pvz, double qvz, double f)
{
	double gp, g0, gm, ddfg;
	gp = gamfqf(typ, pvz, qvz, f, DeltaF);
	g0 = gamfqf(typ, pvz, qvz, f, 0);
	gm = gamfqf(typ, pvz, qvz, f, -DeltaF);
	ddfg = (gp-2.0*g0+gm)/DeltaF/DeltaF;
	return ddfg;
}

int afdisc(Fstype typ, char pfqf, double pvz, double qvz, 
           double *afvec)
{
	double mindf, minda;
	double g, g0;
	double a, a0, da, da0;
	double f, f0, df, df0;
	int max = 2;
	int numwv;
/*
	fprintf(stderr, "afdisc(%c, %c, %g, %g, ...) :\n", typ,pfqf,pvz,qvz); 
*/

	numwv = 0;
	max = Vlsg;
	minda = PI/2.0/((double) max);
	mindf = MaxWvpA/((double) max);
	switch(typ)
	{
	case FS_A: case FS_B: case FS_C: case FS_D: case FS_E:
		a0 = AlphaStart;
		g0 = -0.5*ddagamalph(typ, pvz, qvz, a0);
		da0 = sqrt(GFDiff/g0);
		if(da0 < minda) da0 = minda;
/*
		a = a0;
*/
		a = a0-da0/2.0;
		while(a>=minda/2.0)
		{
			afvec[numwv] = a; ++numwv;
			g = -0.5*ddagamalph(typ, pvz, qvz, a);
			da = sqrt(GFDiff/g);
			if(da < minda) da = minda;
			a = a-da;
		}
/*
		a = a0+da0;
*/
		a = a0+da0/2.0;
		while(a<=PI/2.0-minda/2.0)
		{
			afvec[numwv] = a; ++numwv;
			g = -0.5*ddagamalph(typ, pvz, qvz, a);
			da = sqrt(GFDiff/g);
			if(da < minda) da = minda;
			a = a+da;
		}
		return numwv;

	case FS_F: case FS_G: case FS_H: case FS_I:
		switch(pfqf)
		{
		case 'p': 
			f0 = FStart;
			g0 = -0.5*ddfgamfpf(typ, pvz, qvz, f0);
			df0 = sqrt(GFDiff/g0);
			if(df0 < mindf) df0 = mindf;
/*
			f = f0;
*/
			f = f0-df0/2.0;
			while(f>=mindf/2.0)
			{
				afvec[numwv] = f; ++numwv;
				g = -0.5*ddfgamfpf(typ, pvz, qvz, f);
				df = sqrt(GFDiff/g);
				if(df < mindf) df = mindf;
				f = f-df;
			}
/*
			f = f0+df0;
*/
			f = f0+df0/2.0;
			while(f<=MaxWvpA-mindf/2.0)
			{
				afvec[numwv] = f; ++numwv;
				g = -0.5*ddfgamfpf(typ, pvz, qvz, f);
				df = sqrt(GFDiff/g);
				if(df < mindf) df = mindf;
				f = f+df;
			}
			return numwv;
		case 'q': 
			f0 = FStart;
			g0 = -0.5*ddfgamfqf(typ, pvz, qvz, f0);
			df0 = sqrt(GFDiff/g0);
			if(df0 < mindf) df0 = mindf;
/*
			f = f0;
*/
			f = f0-df0/2.0;
			while(f>=mindf/2.0)
			{
				afvec[numwv] = f; ++numwv;
				g = -0.5*ddfgamfqf(typ, pvz, qvz, f);
				df = sqrt(GFDiff/g);
				if(df < mindf) df = mindf;
				f = f-df;
			}
/*
			f = f0+df0;
*/
			f = f0+df0/2.0;
			while(f<=MaxWvpA-mindf/2.0)
			{
				afvec[numwv] = f; ++numwv;
				g = -0.5*ddfgamfqf(typ, pvz, qvz, f);
				df = sqrt(GFDiff/g);
				if(df < mindf) df = mindf;
				f = f+df;
			}
			return numwv;
		}
	}
	return 0;
}