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

/*
 * wmmdisc.cpp
 * spectral discretisation procedures --- set up the trial function set
 */

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"complex.h"
#include"matrix.h"
#include"wmmmat.h"
#include"waveg.h"
#include"maxweq.h"
#include"wmmtrif.h"
#include"mode.h"
#include"wmmfld.h"
#include"wmmarr.h"
#include"wmm.h"
#include"wmmgam.h"
#include"wmmpol.h"
#include"wmmglob.h"
#include"wmmdisc.h"

/* error message */
void wmmdiscerror(const char *msg)
{
	fprintf(stderr, "WMMDISC: ERROR: %s\n", msg);
	exit(1);
}

/* ................................................................... */

/* set trial function parameters, depending on propagation constant */ 

void settfpars()
{
	int mcp, cp, l, m, j;
	double e;
	Rect r;
/*
fprintf(stderr, "\n settfpars: ");
*/
	mcp = 0; if(Pol == VEC) mcp = 1;
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
			e = Wg.eps(l,m);
			r = Wg.rectbounds(l,m);
			for(cp=0; cp<=mcp; ++cp)
			{
				for(j=0; j<=Phi.ntf(cp,l,m)-1; ++j)
					(*Phi(cp,l,m,j)).initdeppars(
					          Beta, e, K0, Pol, r);
			}
		}
	}
	return;
}

/* ................................................................... */

/* reduce spectral discretization, enforce numerically positive
   normalization matrices */
void reducespecdisc(double b)
{
	int l, m, j, dim, offs, eq;
	int maxi;
	int maxnaf = 0;
	double maxa;
	double *snmat;
	double *ev;

	int r;
	int mn;

	Beta = b;
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
			if(Pol != VEC) mn = Phi.ntf(0,l,m);
			else           mn = Phi.ntf(0,l,m)+Phi.ntf(1,l,m);
			if(mn > maxnaf) maxnaf = mn;
		}
	}
	settfpars();
	snmat = new double[maxnaf*maxnaf];
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
			if(Swf(l,m) != 0.0)
			{
				matto0(snmat, maxnaf, maxnaf);
				dim = fillsnmat(l,m, snmat, 0);
				r = choldc(snmat, dim, 0);
				while(r != 0)
				{
					matto0(snmat, maxnaf, maxnaf);
					dim = fillsnmat(l, m, snmat, 0);
					ev = new double[dim];
					smallesteigval(snmat, dim, 1, ev, 0, 0);
					maxa = fabs(ev[0]);
					maxi = 0;
					for(j=1; j<=dim-1; ++j)
					{
						if(fabs(ev[j]) > maxa)
						{
							maxa = fabs(ev[j]);
							maxi = j;
						}
					}
					delete[] ev;
					eq = 0; j = 0; 
					offs = (*Phi(0,l,m,0)).Amsi;
					while(eq == 0 && j<=Phi.ntf(0,l,m)-1)
					{
						if((*Phi(0,l,m,j)).Amsi-offs 
							== maxi) 
							eq = 1;
						else ++j;
					}
					if(eq == 1)
					{
					   Phi.removetrifun(
						(*Phi(0,l,m,j)).Amsi );
					   --AmDim;
					}
					else
					{
					   if(Pol != VEC)
					      wmmdiscerror("reducepqdisc 1");
					   eq = 0; j = 0; 
					   while(eq == 0 && j<=Phi.ntf(1,l,m)-1)
					   {
					   	if((*Phi(1,l,m,j)).Amsi-offs 
							== maxi) 
							eq = 1;
						else ++j;
					   }
					   if(eq != 1)
					      wmmdiscerror("reducepqdisc 2");
					   Phi.removetrifun(
						(*Phi(1,l,m,j)).Amsi );
					   --AmDim;
					}
					matto0(snmat, maxnaf, maxnaf);
					dim = fillsnmat(l, m, snmat, 0);
					r = choldc(snmat, dim, 0);
				}
			}
					
		}
	}
	delete[] snmat;
	return;
}

/* ................................................................... */

/* the spectral discretization procedure */

/* maximum number of trial functions per rectangle */
#define MNTFPR 1000

/* helper function */
int topqdc(int npq, Fstype typ, char pfqf, int pv, int qv, double *afv,
	    Trifun* tf, int ntf)
{
	int j;
	Trifun *fun;
	if(ntf+npq >= MNTFPR) 
		wmmdiscerror("MNTFPR overflow !!!"); 
	for(j=0; j<=npq-1; ++j)
	{
	   fun = &tf[ntf+j];
           (*fun).Amsi   = AmDim;		
	   (*fun).AFp    = afv[j];
	   (*fun).Pvz    = pv;
	   (*fun).Qvz    = qv;
	   (*fun).Svz    = '+';
	   (*fun).PfQf   = pfqf;
	   (*fun).Xr     = GPxr;
	   (*fun).Yr     = GPyr;
	   (*fun).GfTyp  = typ;
	   (*fun).Gfcof  = 1.0;
	   (*fun).Amplit = 1.0;
	   ++AmDim;
	}
	ntf += npq;
	return ntf;
}

/* helper function */
int addidgf( Fstype typ, 
	     char pfqf, int pv, int qv, double afp, 
	     int si, double xr, double yr, char sv,
	     Trifun* tf, int ntf)
{
	Trifun *fun;
	if(ntf+1 >= MNTFPR) 
		wmmdiscerror("MNTFPR overflow !!!"); 
	fun = &tf[ntf];
        (*fun).Amsi   = si;		
	(*fun).AFp    = afp;
	(*fun).Pvz    = pv;
	(*fun).Qvz    = qv;
	(*fun).Svz    = sv;
	(*fun).PfQf   = pfqf;
	(*fun).Xr     = xr;
	(*fun).Yr     = yr;
	(*fun).GfTyp  = typ;
	(*fun).Gfcof  = 1.0;
	(*fun).Amplit = 1.0;
	return ntf+1;
}

/* the spectral discretization procedure */
void specdisc(double sdm, double b0, double b1, double b2)
{
	int l, m, mm, j;
	int cp, mcp;
	int v, npq, st;
	Fstype typ;
	char sv;
	double yr;
	double w;
	double afv[MTFPSET];

	Trifun tf[MNTFPR];
	Trifun *fun;
	int ntf, ontf; 
	Rect r;
/*
fprintf(stderr, "specdisc(%g, %.16g, %.16g, %.16g)\n", sdm, b0, b1, b2);
*/
/*
fprintf(stderr, "sdm: %g  db: %.16g\n", sdm, b2-b0);
*/

	GFDiff  =     sdm*Par.fin_d_alpha  +(1.0-sdm)*Par.ini_d_alpha;
	MaxWvpA =     sdm*Par.fin_alpha_max+(1.0-sdm)*Par.ini_alpha_max;
	Vlsg    = int(sdm*Par.fin_N_alpha  +(1.0-sdm)*Par.ini_N_alpha);

	Beta = b1;
	SDBeta = Beta;
	AmDim = 0;
	mcp = 0;
	if(Pol == VEC) mcp = 1;
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
		    for(cp=0; cp<=mcp; ++cp)
		    {
			ntf = 0;
			st = Styp(l,m);
			if(cp == 1)
			{
				switch(st)
				{
				case VSYM: st = VASY; break;
				case VASY: st = VSYM; break;
				case IDSY: st = IDAS; break;
				case IDAS: st = IDSY; break;
				}
			}
if(st == CALC || st == VSYM || st == VASY)
{
			w = K0*K0*Wg.eps(l,m)-Beta*Beta; 
			if(w>=0.0) { v = 1; w = sqrt(w);}
			else { v = -1; w = sqrt(-w);}
			r = Wg.rectbounds(l,m);	
			GPxa = r.x0;
			GPxb = r.x1;
			GPya = r.y0;
			GPyb = r.y1;
			GPw  = w;
			switch(Rtyp(l,m))
			{
			case 0: if(v == 1)
				{
				  GPxr = (GPxa+GPxb)/2.0;
				  GPyr = (GPya+GPyb)/2.0;
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_A, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_A, '-', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_B, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_B, '-', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_C, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_C, '-', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_D, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_D, '-', 1, 1, afv,
					tf, ntf);
			    }

				  GPyr = (GPya+GPyb)/2.0; 
				  GPxr = GPxb;
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_F, 'p', 1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'p', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_G, 'p', 1,-1, afv);
                                  ntf = topqdc(npq, FS_G, 'p', 1,-1, afv,
					tf, ntf);
			    }
				  GPxr = GPxa;
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_F, 'p',-1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'p',-1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_G, 'p',-1,-1, afv);
                                  ntf = topqdc(npq, FS_G, 'p',-1,-1, afv,
					tf, ntf);
			    }

				  GPxr = (GPxa+GPxb)/2.0; 
				  GPyr = GPyb;
			    if(st == CALC || st == VASY || st == VSYM)
			    {
				  npq = afdisc(FS_H, 'q', 1, 1, afv);
                                  ntf = topqdc(npq, FS_H, 'q', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY || st == VSYM)
			    {
				  npq = afdisc(FS_I, 'q',-1, 1, afv);
                                  ntf = topqdc(npq, FS_I, 'q',-1, 1, afv,
					tf, ntf);
			    }
				  GPyr = GPya;
			    if(st == CALC)
			    {
				  npq = afdisc(FS_H, 'q', 1,-1, afv);
                                  ntf = topqdc(npq, FS_H, 'q', 1,-1, afv,
					tf, ntf);
			    }
			    if(st == CALC)
			    {
				  npq = afdisc(FS_I, 'q',-1,-1, afv);
                                  ntf = topqdc(npq, FS_I, 'q',-1,-1, afv,
					tf, ntf);
			    }
				}
				else	
				{
				  GPxr = GPxb; GPyr = GPyb;
			    if(st == CALC || st == VSYM || st == VASY)
			    {
				  npq = afdisc(FS_E, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_E, '-', 1, 1, afv,
					tf, ntf);
			    }
				  GPxr = GPxb; GPyr = GPya;
			    if(st == CALC)
			    {
				  npq = afdisc(FS_E, '-', 1,-1, afv);
                                  ntf = topqdc(npq, FS_E, '-', 1,-1, afv,
					tf, ntf);
			    }
				  GPxr = GPxa; GPyr = GPyb;
			    if(st == CALC || st == VSYM || st == VASY)
			    {
				  npq = afdisc(FS_E, '-',-1, 1, afv);
                                  ntf = topqdc(npq, FS_E, '-',-1, 1, afv,
					tf, ntf);
			    }
				  GPxr = GPxa; GPyr = GPya;
			    if(st == CALC)
			    {
				  npq = afdisc(FS_E, '-',-1,-1, afv);
                                  ntf = topqdc(npq, FS_E, '-',-1,-1, afv,
					tf, ntf);
			    }

				  GPxr = (GPxa+GPxb)/2.0; 
				  GPyr = GPyb;
			    if(st == CALC || st == VSYM || st == VASY)
			    {
				  npq = afdisc(FS_H, 'p', 1, 1, afv);
                                  ntf = topqdc(npq, FS_H, 'p', 1, 1, afv,
					tf, ntf);
			    }
				  GPyr = GPya;
			    if(st == CALC)
			    {
				  npq = afdisc(FS_H, 'p', 1,-1, afv);
                                  ntf = topqdc(npq, FS_H, 'p', 1,-1, afv,
					tf, ntf);
			    }
				  GPyr = GPyb;
			    if(st == CALC || st == VSYM || st == VASY)
			    {
				  npq = afdisc(FS_I, 'p',-1, 1, afv);
                                  ntf = topqdc(npq, FS_I, 'p',-1, 1, afv,
					tf, ntf);
			    }
				  GPyr = GPya;
			    if(st == CALC)
			    {
				  npq = afdisc(FS_I, 'p',-1,-1, afv);
                                  ntf = topqdc(npq, FS_I, 'p',-1,-1, afv,
					tf, ntf);
			    }

				  GPyr = (GPya+GPyb)/2.0; 
				  GPxr = GPxb;
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_F, 'q', 1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'q', 1, 1, afv,
					tf, ntf);
			    }
				  GPxr = GPxa;
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_F, 'q',-1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'q',-1, 1, afv,
					tf, ntf);
			    }
				  GPxr = GPxb;
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_G, 'q', 1,-1, afv);
                                  ntf = topqdc(npq, FS_G, 'q', 1,-1, afv,
					tf, ntf);
			    }
				  GPxr = GPxa;
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_G, 'q',-1,-1, afv);
                                  ntf = topqdc(npq, FS_G, 'q',-1,-1, afv,
					tf, ntf);
			    }
				}
				break;
			case 1: if(v==1)
				{
				  GPyr = (GPya+GPyb)/2.0; 
				  GPxr = GPxa;
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_F, 'p',-1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'p',-1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_G, 'p',-1, 1, afv);
                                  ntf = topqdc(npq, FS_G, 'p',-1, 1, afv,
					tf, ntf);
			    }
				}
				else	
				{
				  GPxr = GPxa; GPyr = GPyb;
			    if(st == CALC || st == VSYM || st == VASY) 
			    {
				  npq = afdisc(FS_E, '-',-1, 1, afv);
                                  ntf = topqdc(npq, FS_E, '-',-1, 1, afv,
					tf, ntf);
			    }
				  GPxr = GPxa; GPyr = GPya;
			    if(st == CALC)
			    {
				  npq = afdisc(FS_E, '-',-1,-1, afv);
                                  ntf = topqdc(npq, FS_E, '-',-1,-1, afv,
					tf, ntf);
			    }

				  GPyr = (GPya+GPyb)/2.0; 
				  GPxr = GPxa; 
			    if(st == CALC || st == VSYM)
			    {
				  npq = afdisc(FS_F, 'q',-1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'q',-1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY)
			    {
				  npq = afdisc(FS_G, 'q',-1, 1, afv);
                                  ntf = topqdc(npq, FS_G, 'q',-1, 1, afv,
					tf, ntf);
			    }
				}
				break;
			case 2: 
				GPxr = GPxa; GPyr = GPya;
				npq = afdisc(FS_E, '-',-1,-1, afv);
                                ntf = topqdc(npq, FS_E, '-',-1,-1, afv,
					tf, ntf);
				break;
			case 3: if(v==1)
				{
				  GPxr = (GPxa+GPxb)/2.0; 
				  GPyr = GPya;
				  npq = afdisc(FS_H, 'q', 1,-1, afv);
                                  ntf = topqdc(npq, FS_H, 'q', 1,-1, afv,
					tf, ntf);
				  npq = afdisc(FS_I, 'q', 1,-1, afv);
                                  ntf = topqdc(npq, FS_I, 'q', 1,-1, afv,
					tf, ntf);
				}
				else	
				{
				  GPxr = GPxb; GPyr = GPya;
				  npq = afdisc(FS_E, '-', 1,-1, afv);
                                  ntf = topqdc(npq, FS_E, '-', 1,-1, afv,
					tf, ntf);
				  GPxr = GPxa; GPyr = GPya;
				  npq = afdisc(FS_E, '-',-1,-1, afv);
                                  ntf = topqdc(npq, FS_E, '-',-1,-1, afv,
					tf, ntf);

				  GPxr = (GPxa+GPxb)/2.0; 
				  GPyr = GPya; 
				  npq = afdisc(FS_H, 'p', 1,-1, afv);
                                  ntf = topqdc(npq, FS_H, 'p', 1,-1, afv,
					tf, ntf);
				  npq = afdisc(FS_I, 'p', 1,-1, afv);
                                  ntf = topqdc(npq, FS_I, 'p', 1,-1, afv,
					tf, ntf);
				}
				break;
			case 4: 
				GPxr = GPxb; GPyr = GPya;
				npq = afdisc(FS_E, '-', 1,-1, afv);
                                ntf = topqdc(npq, FS_E, '-', 1,-1, afv,
					tf, ntf);
				break;
			case 5: if(v==1)
				{
				  GPyr = (GPya+GPyb)/2.0; 
				  GPxr = GPxb;
			    if(st == CALC || st == VSYM) 
			    {
				  npq = afdisc(FS_F, 'p', 1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'p', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY) 
			    {
				  npq = afdisc(FS_G, 'p', 1, 1, afv);
                                  ntf = topqdc(npq, FS_G, 'p', 1, 1, afv,
					tf, ntf);
			    }
				}
				else	
				{
				  GPxr = GPxb; GPyr = GPyb;
			    if(st == CALC || st == VSYM || st == VASY) 
			    {
				  npq = afdisc(FS_E, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_E, '-', 1, 1, afv,
					tf, ntf);
			    }
				  GPxr = GPxb; GPyr = GPya;
			    if(st == CALC) 
			    {
				  npq = afdisc(FS_E, '-', 1,-1, afv);
                                  ntf = topqdc(npq, FS_E, '-', 1,-1, afv,
					tf, ntf);
			    }

				  GPyr = (GPya+GPyb)/2.0; 
				  GPxr = GPxb; 
			    if(st == CALC || st == VSYM) 
			    {
				  npq = afdisc(FS_F, 'q', 1, 1, afv);
                                  ntf = topqdc(npq, FS_F, 'q', 1, 1, afv,
					tf, ntf);
			    }
			    if(st == CALC || st == VASY) 
			    {
				  npq = afdisc(FS_G, 'q', 1, 1, afv);
                                  ntf = topqdc(npq, FS_G, 'q', 1, 1, afv,
					tf, ntf);
			    }
				}
				break;
			case 6: 
				GPxr = GPxb; GPyr = GPyb;
				npq = afdisc(FS_E, '-', 1, 1, afv);
                                ntf = topqdc(npq, FS_E, '-', 1, 1, afv,
					tf, ntf);
				break;
			case 7: if(v==1)
				{
				  GPxr = (GPxa+GPxb)/2.0; 
				  GPyr = GPyb;
				  npq = afdisc(FS_H, 'q', 1, 1, afv);
                                  ntf = topqdc(npq, FS_H, 'q', 1, 1, afv,
					tf, ntf);
				  npq = afdisc(FS_I, 'q', 1, 1, afv);
                                  ntf = topqdc(npq, FS_I, 'q', 1, 1, afv,
					tf, ntf);
				}
				else	
				{
				  GPxr = GPxb; GPyr = GPyb;
				  npq = afdisc(FS_E, '-', 1, 1, afv);
                                  ntf = topqdc(npq, FS_E, '-', 1, 1, afv,
					tf, ntf);
				  GPxr = GPxa; GPyr = GPyb;
				  npq = afdisc(FS_E, '-',-1, 1, afv);
                                  ntf = topqdc(npq, FS_E, '-',-1, 1, afv,
					tf, ntf);

				  GPxr = (GPxa+GPxb)/2.0; 
				  GPyr = GPyb; 
				  npq = afdisc(FS_H, 'p', 1, 1, afv);
                                  ntf = topqdc(npq, FS_H, 'p', 1, 1, afv,
					tf, ntf);
				  npq = afdisc(FS_I, 'p', 1, 1, afv);
                                  ntf = topqdc(npq, FS_I, 'p', 1, 1, afv,
					tf, ntf);
				}
				break;
			case 8: 
				GPxr = GPxa; GPyr = GPyb;
				npq = afdisc(FS_E, '-',-1, 1, afv);
                                ntf = topqdc(npq, FS_E, '-',-1, 1, afv,
					tf, ntf);
				break;
			}
}
if(st == VSYM || st == VASY)
{
			if(st == VSYM) sv = '+'; else sv = '-';
			ontf = ntf;
			for(j=0; j<=ontf-1; ++j)
			{
				fun = &tf[j];
				typ = (*fun).GfTyp;
				if(typ == FS_E || typ == FS_H || typ == FS_I)
				{
					ntf = addidgf(typ, (*fun).PfQf, 
						     (*fun).Pvz, -(*fun).Qvz, 
						     (*fun).AFp, (*fun).Amsi, 
						     (*fun).Xr, r.y0, sv,
						     tf, ntf);
				}
			}
}
if(st == VSYM || st == VASY || st == CALC)
{
			Phi.insert(tf, ntf, cp, l, m);
}
		     }
		}
	}

	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
		    for(cp=0; cp<=mcp; ++cp)
		    {
			st = Styp(l,m);
			if(cp == 1)
			{
				switch(st)
				{
				case VSYM: st = VASY; break;
				case VASY: st = VSYM; break;
				case IDSY: st = IDAS; break;
				case IDAS: st = IDSY; break;
				}
			}
if(st == IDSY || st == IDAS)
{
			r = Wg.rectbounds(l, m);
			ntf = 0;
			if(st == IDSY) sv = '+';
			else sv = '-';
			mm = Wg.ny+1-m; 	
			npq = Phi.ntf(cp, l, mm);
			for(j=0; j<=npq-1; ++j)
			{
				fun = Phi(cp, l, mm, j);
				typ = (*fun).GfTyp;
				if(typ == FS_E || typ == FS_H || typ == FS_I)
				{
					if(-(*fun).Qvz > 0) yr = r.y1;
					else yr = r.y0;
				}
				else yr = 0.5*(r.y0+r.y1);
				ntf = addidgf(typ, (*fun).PfQf, 
					     (*fun).Pvz, -(*fun).Qvz, 
					     (*fun).AFp, (*fun).Amsi, 
					     (*fun).Xr, yr, sv, 
					     tf, ntf);
			}
			Phi.insert(tf, ntf, cp, l, m);
}
		    }
		}
	}
	settfpars();	

	j = AmDim;
	fprintf(stderr, " # %d ", AmDim);
	if(b0 == b2) reducespecdisc(b1);
	else
	{
		reducespecdisc(b1);
		reducespecdisc(b0);
		reducespecdisc(b2);
	}
	if(AmDim < j) fprintf(stderr, "> %d\n", AmDim);
	else          fprintf(stderr, "\n");
	return;
}

/* save spectral discretization for diagnosis purposes */
/* helper function: Fstype -> char */
char fstoch(Fstype t) { return 'A'+char(int(t)); }
void diagspecdisc()
{
	FILE *dat;
	int l, m, j;
	Trifun *fun;
	int cp, mcp;
	mcp = 0;
	if(Pol == VEC) mcp = 1;
	fprintf(stderr, ">> pjqj\n");
	dat = fopen("pjqj", "w+");
	fprintf(dat, "beta: %.15g\n", Beta);
	fprintf(dat,"# trial functions:\n");
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
			if(Pol == VEC)
				fprintf(dat, "  [%d,%d] - 0: %d  1: %d\n", 
				l, m, Phi.ntf(0,l,m), Phi.ntf(1,l,m)); 
			else
				fprintf(dat, "  [%d,%d] - %d\n", 
				l, m, Phi.ntf(0,l,m)); 
		}
	}
	fprintf(dat,"amplitude matrix:\n");
	fprintf(dat,"  dimension: %d \n", AmDim);
	fprintf(dat,"trial function parameters:\n");
	for(m=0; m<=Wg.ny+1; ++m)
	{
		for(l=0; l<=Wg.nx+1; ++l)
		{
		    for(cp=0; cp<=mcp; ++cp)
		    {
			fprintf(dat, "[%d,%d] --- cp: %d\n", 
					l, m, cp);
			fprintf(dat, "       #GF: %d  STYP: %d\n", 
				Phi.ntf(cp,l,m), Styp(l,m));
			for(j=0; j<=Phi.ntf(cp,l,m)-1; ++j)
			{
				fun = Phi(cp,l,m,j);
				fprintf(dat, "(%d)(ai: %d)  p: %g q: %g", 
					j, (*fun).Amsi, 
					(*fun).Pj, (*fun).Qj);
				fprintf(dat, 
				   " Gtyp: %c*%g ", 
				   fstoch((*fun).GfTyp), (*fun).Gfcof);
				fprintf(dat, "#: (%g, %g)\n", 
					(*fun).Xr, (*fun).Yr);
			}
			fprintf(dat, "\n");
		    }
		}
	}
	fclose(dat);
	return;
}