/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/*
 * bundle.cpp
 * vectorial QUEP, wave bundles
 */

#include<stdio.h>
#include<stdlib.h>
#include<cmath>
#include"inout.h"
#include"complex.h"
#include"matrix.h"
#include"structure.h"
#include"gengwed.h"
#include"slamode.h"
#include"slamarr.h"
#include"slams.h"
#include"bepfld.h"    
#include"matlvis.h"    
#include"quepfld.h"    
#include"bundle.h"    

/* error message */
void bundleerror(const char *m) 
{
	fprintf(stderr, "\nbundle: %s.\n", m);
	exit(1);
}

/* ------------------------------------------------------------------- 
 * bundles of vQUEP solutions
 */

/* full setup:
   rc: the structure under investigation,
   incoming field: 
      guided mode of polarization p, order o, 
      coming in at border bdr, with unit amplitude, 
      at primary angle of incidence t (^o), 
      beam of full width w (1/e, field),
      represented by na angles between t-dt and t+dt,
      dt determined such that spectral contributions with relative 
      amplitude below dtlev are neglected 
*/
/* vQUEP computations:
      wx: vertical computational window,
      Mx: number of expansion terms per vertical slice,
      wz: horizontal computational window,
      Mx: number of expansion terms per horizontal layer */
Bundle::Bundle(SegWgStruct rc, 
               CBorder bdr, Polarization p, int o, 
               double t, double w, double fcy, double fcz, int na, double dtlev,
               Interval wx, int Mx, Interval wz, int Mz)
{
	CM = 'Q';
	str = rc;
	if(bdr != LEFT) 
		bundleerror("Bundle::Bundle(): bdr != LEFT not yet implemented");
	entryQ = bdr;
	pol = p;
	ord = o;
	theta0 = t;
	bhw = w/2.0;
	y0 = fcy;
	z0 = fcz;

	Nky = na;
	Wky = 2.0/bhw;
	Dky = Wky*sqrt(-log(dtlev));

	cwx = wx;
	nmx = Mx;
	cwz = wz;
	nmz = Mz;

	k0 = val_k0(str(0).lambda);
	Waveguide ewg = str(0);		
	ModeArray ma;
        modeanalysis(ewg, pol, ma, 1);
	if(ma.num < ord+1)
		bundleerror("Bundle::Bundle(): mode order not supported");
	Nin = ma(ord).neff;
	ky0 = k0*Nin*sin(theta0/180.0*PI);
	// m = ma(ord);

	Circuit cstr = str.circuit();

fprintf(stderr, "\n------------- vQUEP-bundle ----------------------------------------- '14 ---\n");
fprintf(stderr, "in: T%c%d  N = %g  ", polCHR(pol), ord, Nin);
fprintf(stderr, "y0 = %g mum  ", y0);
fprintf(stderr, "z0 = %g mum\n", z0);
fprintf(stderr, "width_y = %g mum  ", 2.0*bhw);
fprintf(stderr, "width_cr = %g mum\n", 2.0*bhw*cos(theta0*PI/180.0));
fprintf(stderr, "theta0 = %g^o  ", theta0);
fprintf(stderr, "ky0 = %g /mum\n", ky0);
fprintf(stderr, "ky /mum: [%g :: %g : %g : %g :: %g]_%d\n", 
	ky0-Dky, ky0-Wky, ky0, ky0+Wky, ky0+Dky, Nky);
if(fabs((ky0-Dky)/k0/Nin)<1 && fabs((ky0+Dky)/k0/Nin)<1) 
{
	fprintf(stderr, "theta /^o: [%g :: %g : %g : %g :: %g]_%d\n", 
	asin((ky0-Dky)/k0/Nin)/PI*180.0, 
	asin((ky0-Wky)/k0/Nin)/PI*180.0, 
	asin((ky0    )/k0/Nin)/PI*180.0, 
	asin((ky0+Wky)/k0/Nin)/PI*180.0, 
	asin((ky0+Dky)/k0/Nin)/PI*180.0, Nky);
}
else
{
	fprintf(stderr, "ky_max = %g /mum\n", k0*Nin);
	bundleerror("Bundle::Bundle(): unphysical spectral range");
}
fprintf(stderr, "lambda =  %.10g mum  k0 = %g /mum", str(0).lambda, k0);
if(str.special()) fprintf(stderr, "  (!)\n");
else              fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d,  ", wx.x0, wx.x1, Mx);
fprintf(stderr, "z in (%.10g, %.10g), Mz: %d\n", wz.x0, wz.x1, Mz);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", cstr.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=cstr.nx; ++j) fprintf(stderr, "%6.10g  ", cstr.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    Nz: %d\n", cstr.nz);
fprintf(stderr, "    hz: ");
for(int j=0; j<=cstr.nz; ++j) fprintf(stderr, "%6.10g  ", cstr.hz(j));
fprintf(stderr, "\n");
if(cstr.special) fprintf(stderr, "   eps:\n");
else           fprintf(stderr, "     n:\n");
for(int i=cstr.nx+1; i>=0; --i) 
{
fprintf(stderr,"layer %d:  ", i);
for(int j=0; j<=cstr.nz+1; ++j) fprintf(stderr, "%6.10g  ", cstr.n(i, j));
fprintf(stderr,"\n");
}
fprintf(stderr, "----------------------------------------------------------------------------\n");
	VA.clear();
	FA.clear();
	PA.clear();
	Pin = 0.0;
}
Bundle::Bundle(Circuit rc, 
               CBorder bdr, Polarization p, int o, 
               double t, double w, double fcy, double fcz, int na, double dtlev,
               Interval wx, int Mx, Interval wz, int Mz)
{
	(*this) = Bundle(SegWgStruct(rc), 
	                 bdr, p, o, 
	                 t, w, fcy, fcz, na, dtlev,
	                 wx, Mx, wz, Mz);
}

/* vBEP computations:
      wx: vertical computational window,
      Mx: number of expansion terms per vertical slice */
Bundle::Bundle(SegWgStruct rc, 
               SBorder bdr, Polarization p, int o, 
               double t, double w, double fcy, double fcz, int na, double dtlev,
               Interval wx, int Mx)
{
	CM = 'B';
	str = rc;
	if(bdr != SBLEFT) 
		bundleerror("Bundle::Bundle(): bdr != SBLEFT not yet implemented");
	entryB = bdr;
	pol = p;
	ord = o;
	theta0 = t;
	bhw = w/2.0;
	y0 = fcy;
	z0 = fcz;

	Nky = na;
	Wky = 2.0/bhw;
	Dky = Wky*sqrt(-log(dtlev));

	cwx = wx;
	nmx = Mx;

	k0 = val_k0(str(0).lambda);
	Waveguide ewg = str(0);		
	ModeArray ma;
        modeanalysis(ewg, pol, ma, 1);
	if(ma.num < ord+1)
		bundleerror("Bundle::Bundle(): mode order not supported");
	Nin = ma(ord).neff;
	ky0 = k0*Nin*sin(theta0/180.0*PI);
	// m = ma(ord);

	Circuit cstr = str.circuit();

fprintf(stderr, "\n------------- vBEP-bundle ------------------------------------------ '22 ---\n");
fprintf(stderr, "in: T%c%d  N = %g  ", polCHR(pol), ord, Nin);
fprintf(stderr, "y0 = %g mum  ", y0);
fprintf(stderr, "z0 = %g mum\n", z0);
fprintf(stderr, "width_y = %g mum  ", 2.0*bhw);
fprintf(stderr, "width_cr = %g mum\n", 2.0*bhw*cos(theta0*PI/180.0));
fprintf(stderr, "theta0 = %g^o  ", theta0);
fprintf(stderr, "ky0 = %g /mum\n", ky0);
fprintf(stderr, "ky /mum: [%g :: %g : %g : %g :: %g]_%d\n", 
	ky0-Dky, ky0-Wky, ky0, ky0+Wky, ky0+Dky, Nky);
if(fabs((ky0-Dky)/k0/Nin)<1 && fabs((ky0+Dky)/k0/Nin)<1) 
{
	fprintf(stderr, "theta /^o: [%g :: %g : %g : %g :: %g]_%d\n", 
	asin((ky0-Dky)/k0/Nin)/PI*180.0, 
	asin((ky0-Wky)/k0/Nin)/PI*180.0, 
	asin((ky0    )/k0/Nin)/PI*180.0, 
	asin((ky0+Wky)/k0/Nin)/PI*180.0, 
	asin((ky0+Dky)/k0/Nin)/PI*180.0, Nky);
}
else
{
	fprintf(stderr, "ky_max = %g /mum\n", k0*Nin);
	bundleerror("Bundle::Bundle(): unphysical spectral range");
}
fprintf(stderr, "lambda =  %.10g mum  k0 = %g /mum", str(0).lambda, k0);
if(str.special()) fprintf(stderr, "  (!)\n");
else              fprintf(stderr, "\n");
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "x in (%.10g, %.10g), Mx: %d\n", wx.x0, wx.x1, Mx);
fprintf(stderr, "----------------------------------------------------------------------------\n");
fprintf(stderr, "    Nx: %d\n", cstr.nx);
fprintf(stderr, "    hx: ");
for(int j=0; j<=cstr.nx; ++j) fprintf(stderr, "%6.10g  ", cstr.hx(j));
fprintf(stderr, "\n");
fprintf(stderr, "    Nz: %d\n", cstr.nz);
fprintf(stderr, "    hz: ");
for(int j=0; j<=cstr.nz; ++j) fprintf(stderr, "%6.10g  ", cstr.hz(j));
fprintf(stderr, "\n");
if(cstr.special) fprintf(stderr, "   eps:\n");
else           fprintf(stderr, "     n:\n");
for(int i=cstr.nx+1; i>=0; --i) 
{
fprintf(stderr,"layer %d:  ", i);
for(int j=0; j<=cstr.nz+1; ++j) fprintf(stderr, "%6.10g  ", cstr.n(i, j));
fprintf(stderr,"\n");
}
fprintf(stderr, "----------------------------------------------------------------------------\n");
	VA.clear();
	FA.clear();
	PA.clear();
	Pin = 0.0;
}
Bundle::Bundle(Circuit rc, 
               SBorder bdr, Polarization p, int o, 
               double t, double w, double fcy, double fcz, int na, double dtlev,
               Interval wx, int Mx)
{
	(*this) = Bundle(SegWgStruct(rc), 
	                 bdr, p, o, 
	                 t, w, fcy, fcz, na, dtlev,
	                 wx, Mx);
}

/* field for a single lateral spectral component */
Complex Bundle::spfld(Fcomp cp, double x, double y, double z)
{
	if(CM == 'Q') return qfld.field(cp, x, z)*expmi(ky*(y-y0));
	return               bfld.field(cp, x, z)*expmi(ky*(y-y0));
}

/* stepwise vQUEP: fill views with field data,
   vQUEP calcuations with shifting boundaries,  
   numtr, pviol, bshift: arguments for quepsim_s() */
void Bundle::solve(int numtr, double pviol, double bshift)
{
	if(CM != 'Q') bundleerror("solve(Q): wrong compuational model");
	if(Nky <= 6) bundleerror("solve(Q): Nky too small");
	double kys = 2.0*Dky/((double) Nky);
	for(int sci=0; sci<=Nky; ++sci)
	{
		ky = ky0-Dky+((double) sci)*kys;
		double theta = asin(ky/k0/Nin)/PI*180.0;
		double kz = sqrt(k0*k0*Nin*Nin-ky*ky);
		Complex amp = expmi(kz*(cwz.x0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		fprintf(stderr, "\n[%d] theta = %g ^o, amp = %g\n", sci, theta, amp.abs()); 

		qfld = QuepField(str,  LEFT, pol, 0, amp, theta, cwx, nmx, cwz, nmz);
	     	int succ = qfld.quepsim_s(numtr, pviol, bshift);
		fprintf(stderr, "P_out/P_in = %g", 
			(qfld.Pout(LEFT)+qfld.Pout(RIGHT)
			+qfld.Pout(BOTTOM)+qfld.Pout(TOP))/amp.sqabs());
		if(succ == 0) fprintf(stderr, "; probably corrupt field!\n");
		else 	      fprintf(stderr, "; ok.\n");
		double                       icoeff =       1.0*kys;
		if(sci == 0 || sci == Nky  ) icoeff =   3.0/8.0*kys;
		if(sci == 1 || sci == Nky-1) icoeff =   7.0/6.0*kys; 
		if(sci == 2 || sci == Nky-2) icoeff = 23.0/24.0*kys; 

		for(int n=0; n<=VA.num-1; ++n)
		{ 
			double dv = (VA(n).va_end-VA(n).va_beg)
			            /((double) (VA(n).va_N-1));
			double dh = (VA(n).ha_end-VA(n).ha_beg)
			            /((double) (VA(n).ha_N-1));
			for(int i=0; i<=VA(n).va_N-1; ++i)
			{
				double v = VA(n).va_beg+((double) i)*dv;
				for(int j=0; j<=VA(n).ha_N-1; ++j)
				{
					double h = VA(n).ha_beg+((double) j)*dh;
					double x=0.0;
					double y=0.0;
					double z=0.0;
					VA(n).trf(v, h, x, y, z);
					VA(n).ex(i, j) += spfld(EX, x, y, z)*icoeff;
					VA(n).ey(i, j) += spfld(EY, x, y, z)*icoeff;
					VA(n).ez(i, j) += spfld(EZ, x, y, z)*icoeff;
					VA(n).hx(i, j) += spfld(HX, x, y, z)*icoeff;
					VA(n).hy(i, j) += spfld(HY, x, y, z)*icoeff;
					VA(n).hz(i, j) += spfld(HZ, x, y, z)*icoeff;
				}
			}
		}
		for(int n=0; n<=FA.num-1; ++n)
		{  
			double x = FA(n).x;
			double y = FA(n).y;
			double z = FA(n).z;
			FA(n).ex += spfld(EX, x, y, z)*icoeff;
			FA(n).ey += spfld(EY, x, y, z)*icoeff;
			FA(n).ez += spfld(EZ, x, y, z)*icoeff;
			FA(n).hx += spfld(HX, x, y, z)*icoeff;
			FA(n).hy += spfld(HY, x, y, z)*icoeff;
			FA(n).hz += spfld(HZ, x, y, z)*icoeff;
		}
//		double pamp = exp(-2.0*(ky-ky0)*(ky-ky0)/(Wky*Wky));
		Pin += qfld.Pin(entryQ)*icoeff;
		for(int n=0; n<=PA.num-1; ++n)
		{  
			Complex a = qfld.Aout(PA(n).bdrQ, PA(n).pol, PA(n).ord);
			PA(n).Pout += a.sqabs()*icoeff;
		}
		qfld.clear();
	}
	for(int n=0; n<=FA.num-1; ++n) FA(n).eps = str.eps(FA(n).x, FA(n).z);
	return;
}

/* stepwise vBEP: fill views with field data */
void Bundle::solve()
{
	if(CM != 'B') bundleerror("solve(B): wrong compuational model");
	if(Nky <= 6) bundleerror("solve(B): Nky too small");
	double kys = 2.0*Dky/((double) Nky);
	for(int sci=0; sci<=Nky; ++sci)
	{
		ky = ky0-Dky+((double) sci)*kys;
		double theta = asin(ky/k0/Nin)/PI*180.0;
		double kz = sqrt(k0*k0*Nin*Nin-ky*ky);
		Complex amp = expmi(kz*(cwz.x0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		fprintf(stderr, "\n[%d] theta = %g ^o, amp = %g\n", sci, theta, amp.abs()); 

		bfld = BepField(str, entryB, pol, ord, amp, theta, cwx, nmx);
	     	bfld.bepsim();
		fprintf(stderr, "P_out/P_in = %g", 
			(bfld.Pout(SBLEFT)+bfld.Pout(SBRIGHT))/amp.sqabs());
		fprintf(stderr, "; ok.\n");
		double                       icoeff =       1.0*kys;
		if(sci == 0 || sci == Nky  ) icoeff =   3.0/8.0*kys;
		if(sci == 1 || sci == Nky-1) icoeff =   7.0/6.0*kys; 
		if(sci == 2 || sci == Nky-2) icoeff = 23.0/24.0*kys; 

		for(int n=0; n<=VA.num-1; ++n)
		{ 
			double dv = (VA(n).va_end-VA(n).va_beg)
			            /((double) (VA(n).va_N-1));
			double dh = (VA(n).ha_end-VA(n).ha_beg)
			            /((double) (VA(n).ha_N-1));
			for(int i=0; i<=VA(n).va_N-1; ++i)
			{
				double v = VA(n).va_beg+((double) i)*dv;
				for(int j=0; j<=VA(n).ha_N-1; ++j)
				{
					double h = VA(n).ha_beg+((double) j)*dh;
					double x=0.0;
					double y=0.0;
					double z=0.0;
					VA(n).trf(v, h, x, y, z);
					VA(n).ex(i, j) += spfld(EX, x, y, z)*icoeff;
					VA(n).ey(i, j) += spfld(EY, x, y, z)*icoeff;
					VA(n).ez(i, j) += spfld(EZ, x, y, z)*icoeff;
					VA(n).hx(i, j) += spfld(HX, x, y, z)*icoeff;
					VA(n).hy(i, j) += spfld(HY, x, y, z)*icoeff;
					VA(n).hz(i, j) += spfld(HZ, x, y, z)*icoeff;
				}
			}
		}
		for(int n=0; n<=FA.num-1; ++n)
		{  
			double x = FA(n).x;
			double y = FA(n).y;
			double z = FA(n).z;
			FA(n).ex += spfld(EX, x, y, z)*icoeff;
			FA(n).ey += spfld(EY, x, y, z)*icoeff;
			FA(n).ez += spfld(EZ, x, y, z)*icoeff;
			FA(n).hx += spfld(HX, x, y, z)*icoeff;
			FA(n).hy += spfld(HY, x, y, z)*icoeff;
			FA(n).hz += spfld(HZ, x, y, z)*icoeff;
		}
//		double pamp = exp(-2.0*(ky-ky0)*(ky-ky0)/(Wky*Wky));
		Pin += bfld.Pin(entryB)*icoeff;
		for(int n=0; n<=PA.num-1; ++n)
		{  
			Complex a = bfld.Aout(PA(n).bdrB, PA(n).pol, PA(n).ord);
			PA(n).Pout += a.sqabs()*icoeff;
		}
		bfld.clear();
	}
	for(int n=0; n<=FA.num-1; ++n) FA(n).eps = str.eps(FA(n).x, FA(n).z);
	return;
}

/* vBEP only: double pass solver, consider transmission resonances 
   in output to mode of order oo, with polarization op, at boundary obdr, 
   resolved with a stepsize that is reduced by a factor of srf */
void Bundle::solveR(SBorder obdr, Polarization op, int oo, double srf)
{
	if(CM != 'B') bundleerror("solveR: implemented for vBEP only");
	if(Nky <= 6) bundleerror("solve(B): Nky too small");

	double ka = ky0-Dky;
	double kb = ky0+Dky;
	double kys = (kb-ka)/((double) Nky);
	fprintf(stderr, "\n ++E   ky /mum = [%g, %g]_%d\n", ka, kb, Nky); 
	fprintf(stderr, " ++E theta /^o = [%g, %g]\n", asin(ka/k0/Nin)/PI*180.0, asin(kb/k0/Nin)/PI*180.0);

	Dvector knp, anp;
	knp.strip();
	anp.strip();

	Dvector th, am;
	th.strip();
	am.strip();
	int done=0;
	while(done==0)
	{
		done = next(Interval(ka, kb), kys, kys*srf, th, am, 1, 0.01);
		if(done==0)
		{
			ky = th(th.nel-1);
			double theta = asin(ky/k0/Nin)/PI*180.0;

			double kz = sqrt(k0*k0*Nin*Nin-ky*ky);
			Complex amp = expmi(kz*(cwz.x0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
			fprintf(stderr, "\n ++[E] theta = %g^o, amp = %g\n", theta, amp.abs()); 

			bfld = BepField(str, entryB, pol, ord, amp, theta, cwx, nmx);
	     		bfld.bepsim();
			Complex a;
			a = bfld.Aout(obdr, op, oo);
			double aq = a.sqabs();
			fprintf(stderr, "P_out/P_in = %g.\n", (bfld.Pout(SBLEFT)+bfld.Pout(SBRIGHT))/amp.sqabs());
			fprintf(stderr, "\n ++[E]R theta = %g ^o, |A_out|^2 = %g\n", theta, aq); 
			am.append(aq);
		
			knp.append(ky);
			anp.append(aq);
		}
	}
	sort(knp, anp);
toxyf("expky0", knp, anp);
	double amax = anp.max();
	double da = amax/Nky;
	Dvector knpn;
	knpn.strip();
Dvector anpn;
anpn.strip();
	for(int i=0; i<=knp.nel-2; ++i)
	{
		if(fabs(anp(i+1)-anp(i)) > da)
		{
			int nap = floor(fabs(anp(i+1)-anp(i))/da);
			double s = (knp(i+1)-knp(i))/((double)(nap+1));
			for(int j=1; j<=nap; ++j)
			{
				knpn.append(knp(i)+((double)j)*s);
anpn.append(-1.0);
			}
		}
	}
	Dvector knpf(knp.nel+knpn.nel);
Dvector anpf(anp.nel+anpn.nel);
	if(knpn.nel > 0)
	{
		knpf.setsubvector(knp, 0);
		knpf.setsubvector(knpn, knp.nel);
anpf.setsubvector(anp, 0);
anpf.setsubvector(anpn, anp.nel);
	}
	else 
	{
		knpf = knp;
anpf = anp;
	}
sort(knpf, anpf);
toxyf("expky", knpf, anpf);

	knpf.sort();
	int nk = knpf.nel;
	fprintf(stderr, "\n +++   ky /mum = [%g, %g]_%d\n", knpf(0), knpf(nk-1), nk); 
	fprintf(stderr, " +++ theta /^o = [%g, %g]\n", asin(knpf(0)/k0/Nin)/PI*180.0, asin(knpf(nk-1)/k0/Nin)/PI*180.0);
	double A = (1.0/sqrt(2.0*PI*sqrt(PI/2.0)))/sqrt(Wky);
	for(int ik=0; ik<=nk-1; ++ik)
	{
		ky = knpf(ik);
		double theta = asin(ky/k0/Nin)/PI*180.0;
		double kz = sqrt(k0*k0*Nin*Nin-ky*ky);
		Complex amp = A*expmi(kz*(cwz.x0-z0))*exp(-(ky-ky0)*(ky-ky0)/(Wky*Wky));
		fprintf(stderr, "\n ++[%d] theta = %g ^o, amp = %g\n", ik, theta, amp.abs()); 

		bfld = BepField(str, entryB, pol, ord, amp, theta, cwx, nmx);
     		bfld.bepsim();
		fprintf(stderr, "P_out/P_in = %g.\n", (bfld.Pout(SBLEFT)+bfld.Pout(SBRIGHT))/amp.sqabs());

		double icoeff = 0.0;
		if(ik == 0)              icoeff = 0.5*(knpf(1)-knpf(0));
		if(ik >=1 && ik <= nk-2) icoeff = 0.5*(knpf(ik+1)-knpf(ik-1));
		if(ik == nk-1)           icoeff = 0.5*(knpf(nk-1)-knpf(nk-2));

		for(int n=0; n<=VA.num-1; ++n)
		{ 
			double dv = (VA(n).va_end-VA(n).va_beg)
			            /((double) (VA(n).va_N-1));
			double dh = (VA(n).ha_end-VA(n).ha_beg)
			            /((double) (VA(n).ha_N-1));
			for(int i=0; i<=VA(n).va_N-1; ++i)
			{
				double v = VA(n).va_beg+((double) i)*dv;
				for(int j=0; j<=VA(n).ha_N-1; ++j)
				{
					double h = VA(n).ha_beg+((double) j)*dh;
					double x=0.0;
					double y=0.0;
					double z=0.0;
					VA(n).trf(v, h, x, y, z);
					VA(n).ex(i, j) += spfld(EX, x, y, z)*icoeff;
					VA(n).ey(i, j) += spfld(EY, x, y, z)*icoeff;
					VA(n).ez(i, j) += spfld(EZ, x, y, z)*icoeff;
					VA(n).hx(i, j) += spfld(HX, x, y, z)*icoeff;
					VA(n).hy(i, j) += spfld(HY, x, y, z)*icoeff;
					VA(n).hz(i, j) += spfld(HZ, x, y, z)*icoeff;
				}
			}
		}
		for(int n=0; n<=FA.num-1; ++n)
		{  
			double x = FA(n).x;
			double y = FA(n).y;
			double z = FA(n).z;
			FA(n).ex += spfld(EX, x, y, z)*icoeff;
			FA(n).ey += spfld(EY, x, y, z)*icoeff;
			FA(n).ez += spfld(EZ, x, y, z)*icoeff;
			FA(n).hx += spfld(HX, x, y, z)*icoeff;
			FA(n).hy += spfld(HY, x, y, z)*icoeff;
			FA(n).hz += spfld(HZ, x, y, z)*icoeff;
		}
		for(int n=0; n<=PA.num-1; ++n)
		{  
			Complex a = bfld.Aout(PA(n).bdrB, PA(n).pol, PA(n).ord);
			PA(n).Pout += a.sqabs()*icoeff;
		}
		Pin += amp.sqabs()*icoeff;
		bfld.clear();
	}
	for(int n=0; n<=FA.num-1; ++n) FA(n).eps = str.eps(FA(n).x, FA(n).z);
	return;
}

/* ------------------------------------------------------------------- 
 * containers for registration of guided modal output power 
 */

/* initialize */
Port::Port()
{
	bdrQ = LEFT;
	bdrB = SBLEFT;
	pol = TE;
	ord = 0;
	id = '?';
	Pout = 0.0;
}
Port::Port(CBorder b, Polarization p, int o, char d)
{
	bdrQ = b;
	if(bdrQ == LEFT)  bdrB = SBLEFT;
	if(bdrQ == RIGHT) bdrB = SBRIGHT;
	pol = p;
	ord = o;
	id  = d;
	Pout = 0.0;
}
Port::Port(SBorder b, Polarization p, int o, char d)
{
	bdrB = b;
	if(bdrB == SBLEFT)  bdrQ = LEFT;
	if(bdrB == SBRIGHT) bdrQ = RIGHT;
	pol = p;
	ord = o;
	id  = d;
	Pout = 0.0;
}
   
/* ------------------------------------------------------------------- 
 * an array of containers for registration of guided modal output power 
 */
   
/* initialize */
PortArray::PortArray()
          : num(0)
{
	pvec = NULL;
}

/* destroy */
PortArray::~PortArray()
{
	if(pvec != NULL) delete[] pvec; 
	pvec = NULL;
	num = 0; 
}

/* copy constructor */
PortArray::PortArray(const PortArray& ma)
	  : num(ma.num)
{
	pvec = new Port[num];
	Port* ap = ma.pvec;
	Port* mp =    pvec;
	for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
}
	
/* assignment */
PortArray& PortArray::operator=(const PortArray& ma)
{
	if(this != &ma)
	{
		if(pvec != NULL) delete[] pvec;
		num = ma.num;
		pvec = new Port[num];
		Port *ap = ma.pvec;
		Port *mp =    pvec;
		for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
	}
	return *this;
}

/* delete all Port entries */
void PortArray::clear()
{
	if(pvec != NULL) delete[] pvec;
	pvec = NULL;
	num = 0; 
}

/* subscripting */
Port& PortArray::operator() (int i)  
{ 
	if(i<0 || i>=num) bundleerror("PortArray: () out of range");
	return pvec[i]; 
}
Port PortArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) bundleerror("PortArray: () out of range");
	return pvec[i]; 
}

/* add a port */
void PortArray::add(Port m)
{
	Port *avec;
	avec = new Port[num+1];
	Port* ap = avec;
	Port* mp = pvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	*ap = m;
	if(pvec != NULL) delete[] pvec;
	pvec = avec;
	++num;
	return;
}

/* ------------------------------------------------------------------- 
 * an array of containers for discretized fields 
 */
   
/* initialize */
ViewArray::ViewArray()
          : num(0)
{
	vvec = NULL;
}

/* destroy */
ViewArray::~ViewArray()
{
	if(vvec != NULL) delete[] vvec; 
	vvec = NULL;
	num = 0; 
}

/* copy constructor */
ViewArray::ViewArray(const ViewArray& ma)
	  : num(ma.num)
{
	vvec = new View[num];
	View* ap = ma.vvec;
	View* mp =    vvec;
	for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
}
	
/* assignment */
ViewArray& ViewArray::operator=(const ViewArray& ma)
{
	if(this != &ma)
	{
		if(vvec != NULL) delete[] vvec;
		num = ma.num;
		vvec = new View[num];
		View *ap = ma.vvec;
		View *mp =    vvec;
		for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
	}
	return *this;
}

/* delete all View entries */
void ViewArray::clear()
{
	if(vvec != NULL) delete[] vvec;
	vvec = NULL;
	num = 0; 
}

/* subscripting */
View& ViewArray::operator() (int i)  
{ 
	if(i<0 || i>=num) bundleerror("ViewArray: () out of range");
	return vvec[i]; 
}
View ViewArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) bundleerror("ViewArray: () out of range");
	return vvec[i]; 
}

/* add a view */
void ViewArray::add(View m)
{
	View *avec;
	avec = new View[num+1];
	View* ap = avec;
	View* mp = vvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	*ap = m;
	if(vvec != NULL) delete[] vvec;
	vvec = avec;
	++num;
	return;
}

/* ------------------------------------------------------------------- 
 * containers for discretized fields, plot routines 
 */

/* initialize */
View::View()
{
	type = XZ;
	ha_beg = -1.0;;
	ha_end = 1.0;
	ha_N   = 10;
	va_beg = -1.0;
	va_end = 1.0;
	va_N   = 10;
	pos = 0.0;
	ex = Cmatrix(va_N, ha_N); ex.init(CC0);
	ey = Cmatrix(va_N, ha_N); ey.init(CC0);
	ez = Cmatrix(va_N, ha_N); ez.init(CC0);
	hx = Cmatrix(va_N, ha_N); hx.init(CC0);
	hy = Cmatrix(va_N, ha_N); hy.init(CC0);
	hz = Cmatrix(va_N, ha_N); hz.init(CC0);
	va_chr = 'x';
	ha_chr = 'z';
	st = SegWgStruct();
	pid0 = '-';
	pid1 = '-';
}

/* initialize: supply discretization information, setup matrices */
View::View(Vtype t, SegWgStruct s, 
           double v_beg, double v_end, int v_N, 
           double h_beg, double h_end, int h_N, double p, 
           char ext0, char ext1)
{
	type = t;
	ha_beg = h_beg;
	ha_end = h_end;
	ha_N   = h_N;
	va_beg = v_beg;
	va_end = v_end;
	va_N   = v_N;
	pos = p;
	ex = Cmatrix(va_N, ha_N); ex.init(CC0);
	ey = Cmatrix(va_N, ha_N); ey.init(CC0);
	ez = Cmatrix(va_N, ha_N); ez.init(CC0);
	hx = Cmatrix(va_N, ha_N); hx.init(CC0);
	hy = Cmatrix(va_N, ha_N); hy.init(CC0);
	hz = Cmatrix(va_N, ha_N); hz.init(CC0);
	va_chr = 'v';
	ha_chr = 'h';
	switch(type)
	{
	case XZ: va_chr = 'x'; ha_chr = 'z'; break;
	case YZ: va_chr = 'y'; ha_chr = 'z'; break;
	case XY: va_chr = 'x'; ha_chr = 'y'; break;
	}
	Waveguide wg;
	double locn;
	switch(type)
	{
	case XZ: 
		st = s;
		break;
	case YZ: 
		st = s;
		for(int j=0; j<=s.nz+1; ++j)
		{
			locn = s(j).n(s(j).layeridx(pos)); 
			wg = Waveguide(1);
			wg.hx(0) = va_beg+(va_end-va_beg)*1.0/3.0; 
			wg.hx(1) = va_beg+(va_end-va_beg)*2.0/3.0; 
			wg.n.init(locn);
			wg.lambda = s(j).lambda;
			st(j) = wg;
		}
		break;
	case XY: 
		st = SegWgStruct(1);
		st.hz(0) = ha_beg+(ha_end-ha_beg)*1.0/3.0; 
		st.hz(1) = ha_beg+(ha_end-ha_beg)*2.0/3.0; 
		wg = s(s.segidx(pos));
		st.init(wg);
		break;
	}
	pid0 = ext0;
	pid1 = ext1;
}

/* transform v, h, pos to x, y, z coordinates */
void View::trf(double v, double h, double& x, double& y, double& z) const
{
	switch(type)
	{
	case XZ: 
		x = v; y = pos; z = h; return;
		break;
	case YZ: 
		x = pos; y = v; z = h; return;
		break;
	case XY: 
		x = v; y = h; z = pos; return;
		break;
	}
	return;
}

/* visualization: output to MATLAB .m files 
 * ext0, ext1: filename id characters for the m.file */

/* image plot corresponding to field component cp
   cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ 
   foa: ORG, MOD, SQR, REP, IMP (ORG = REP) */
void View::plot(Fcomp cp, Afo foa) const 
{
	FILE *dat;
	char name[15] = "__pl_____I.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[4] = afochr(foa);
        name[5] = fldchr(cp);
        name[6] = cpchr(cp);      
	name[7] = pid0;
	name[8] = pid1;
	fprintf(stderr, "plot_%c%c([%g, %g] x [%g, %g] @ %g) >> %s\n", 
	                va_chr, ha_chr, va_beg, va_end, ha_beg, ha_end, pos, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Dmatrix fld;
	switch(cp)
	{
		case EX: fld = realize(ex, foa); break;
		case EY: fld = realize(ey, foa); break;
		case EZ: fld = realize(ez, foa); break;
		case HX: fld = realize(hx, foa); break;
		case HY: fld = realize(hy, foa); break;
		case HZ: fld = realize(hz, foa); break;
		default: bundleerror("View::plot: cp"); break;
	}
	mlout_fld(dat, va_N, ha_N, cp, fld);
	name[10] = 0; 
	char desc[10];
	afocpdesc(foa, cp, desc);
	mlout_genimage(va_chr, ha_chr, cp, foa, name, desc, dat); 
	if(foa == MOD || foa == SQR) mlout_lincolormap(dat);
	else                         mlout_magcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* image plot of the field "strength", id: one of mE, mH, qE, qH, mW */
void View::plot(FSid id) const
{
	if(va_N <= 1 || ha_N <= 1) 
		bundleerror("View::plot, FSid: va_N, ha_N"); 
	FILE *dat;
	char name[15] = "__pl____I.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[4] = mqchr(id);
	name[5] = idchr(id);
	name[6] = pid0;
	name[7] = pid1;
	fprintf(stderr, "plot_%c%c([%g, %g] x [%g, %g] @ %g) >> %s\n", 
	                va_chr, ha_chr, va_beg, va_end, ha_beg, ha_end, pos, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Dmatrix fld(va_N, ha_N);
	double eps;
	double dv = (va_end-va_beg)/((double) (va_N-1));
	double dh = (ha_end-ha_beg)/((double) (ha_N-1));
	double vp, hp;

	hp = ha_beg;
	for(int j=0; j<=ha_N-1; ++j)
	{
		vp = va_beg;
		for(int i=0; i<=va_N-1; ++i) 
		{
			switch(id)
			{
			case mE:
				fld(i, j) = sqrt( ex(i, j).sqabs()
					         +ey(i, j).sqabs()
					         +ez(i, j).sqabs());
				break;
			case mH:
				fld(i, j) = sqrt( hx(i, j).sqabs()
					         +hy(i, j).sqabs()
					         +hz(i, j).sqabs());
				break;
			case qE:
				fld(i, j) = ( ex(i, j).sqabs()
					     +ey(i, j).sqabs()
					     +ez(i, j).sqabs());
				break;
			case qH:
				fld(i, j) = ( hx(i, j).sqabs()
					     +hy(i, j).sqabs()
					     +hz(i, j).sqabs());
				break;
			case mW:
				eps = st.eps(vp, hp);
				fld(i, j) = 0.25*(    MU0*(hx(i, j).sqabs()
							  +hy(i, j).sqabs()
							  +hz(i, j).sqabs())
						 +EP0*eps*(ex(i, j).sqabs()
							  +ey(i, j).sqabs()
							  +ez(i, j).sqabs()));
				break;
			default:
				bundleerror("plot, FSid?");
				break;
			}
			vp += dv;
		}
		hp += dh;
	}
	mlout_fld(dat, va_N, ha_N, SZ, fld);
	name[9] = 0; 
	char desc[10];
	iddesc(id, desc);
	mlout_genimage(va_chr, ha_chr, SZ, MOD, name, desc, dat); 
	mlout_lowlevcolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* electromagnetic energy density image */
void View::plot() const
{
	if(va_N <= 1 || ha_N <= 1) 
		bundleerror("View::plot, energy density: va_N, ha_N"); 
	FILE *dat;
	char name[15] = "__plW__I.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[5] = pid0;
	name[6] = pid1;
	fprintf(stderr, "plot_%c%c([%g, %g] x [%g, %g] @ %g) >> %s\n", 
	                va_chr, ha_chr, va_beg, va_end, ha_beg, ha_end, pos, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference field");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Dmatrix fld(va_N, ha_N);
	double eps;
	double dv = (va_end-va_beg)/((double) (va_N-1));
	double dh = (ha_end-ha_beg)/((double) (ha_N-1));
	double vp, hp;

	hp = ha_beg;
	for(int j=0; j<=ha_N-1; ++j)
	{
		vp = va_beg;
		for(int i=0; i<=va_N-1; ++i) 
		{
			
			eps = st.eps(vp, hp);
			fld(i, j) = 0.25*(    MU0*(hx(i, j).sqabs()
			                          +hy(i, j).sqabs()
			                          +hz(i, j).sqabs())
                                         +EP0*eps*(ex(i, j).sqabs()
			                          +ey(i, j).sqabs()
			                          +ez(i, j).sqabs()));
			vp += dv;
		}
		hp += dh;
	}
	mlout_fld(dat, va_N, ha_N, SZ, fld);
	name[8] = 0; 
	char desc[3];
	desc[0] = 'W'; desc[1] = 0;
	mlout_genimage(va_chr, ha_chr, SZ, MOD, name, desc, dat); 
	mlout_lincolormap(dat);
	mlout_print(dat, name, 'e'); 
	fclose(dat);
	return;
}

/* animation of the light propagation
   the frames show images of component cp
   at ntfr equidistant times over one time period */
void View::movie(Fcomp cp, int ntfr) const
{
	FILE *dat;
	char name[15] = "__an____M.m";
	double tmax, dt;
        name[0] = va_chr;
        name[1] = ha_chr;
        name[4] = fldchr(cp);
        name[5] = cpchr(cp);      
	name[6] = pid0;
	name[7] = pid1;
	tmax = val_period_T(st(0).lambda);
	dt = tmax/ntfr;
	fprintf(stderr, 
	   "movie_%c%c(%d x [%g, %g] x [%g, %g] @ %g), T=%g fs, dt=%g fs >> %s\n", 
	   va_chr, ha_chr, ntfr, va_beg, va_end, ha_beg, ha_end, pos, tmax, dt, name);
	dat = fopen(name, "w+");
	mlout_title(dat, name, "Interference animation");
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);
	Cmatrix fld;
	switch(cp)
	{
		case EX: fld = ex; break;
		case EY: fld = ey; break;
		case EZ: fld = ez; break;
		case HX: fld = hx; break;
		case HY: fld = hy; break;
		case HZ: fld = hz; break;
		default: bundleerror("View::movie: cp"); break;
	}
	double famp = fld.max();
	mlout_fld(dat, va_N, ha_N, cp, fld.re());
	mlout_fldtore(dat, cp);
	mlout_fld(dat, va_N, ha_N, cp, fld.im());
	mlout_fldtoim(dat, cp);
	name[9] = 0; 
	char desc[10];
	afocpdesc(ORG, cp, desc);
	mlout_genmovie(va_chr, ha_chr, cp, name, desc, dat, ntfr, dt, st(0).lambda, famp); 
	fclose(dat);
	return;
}

/* export full field data ("everything") into a viewer MATLAB file */
void View::viewer() const
{
	FILE *dat;
	char name[13] = "__fld__A.m";
	name[0] = va_chr;
	name[1] = ha_chr;
	name[5] = pid0;
	name[6] = pid1;
	fprintf(stderr, 
	   "viewer_%c%c([%g (%d) %g] x [%g (%d) %g] @ %g) >> %s\n", 
	   va_chr, ha_chr, va_beg, va_N, va_end, ha_beg, ha_N, ha_end, pos, name);
	dat = fopen(name, "w+");
	name[8] = 0; 
	mlout_viewertop(dat, name, va_chr, ha_chr, TE, st(0).lambda, 1);
	mlout_gengeoxz(dat, st, va_beg, va_end, ha_beg, ha_end);
	mlout_mesh(dat, va_chr, va_beg, va_end, va_N, ha_chr, ha_beg, ha_end, ha_N);

	mlout_fld(dat, va_N, ha_N, EX, ex.re()); mlout_fldtore(dat, EX);
	mlout_fld(dat, va_N, ha_N, EX, ex.im()); mlout_fldtoim(dat, EX);
	mlout_fld(dat, va_N, ha_N, EY, ey.re()); mlout_fldtore(dat, EY);
	mlout_fld(dat, va_N, ha_N, EY, ey.im()); mlout_fldtoim(dat, EY);
	mlout_fld(dat, va_N, ha_N, EZ, ez.re()); mlout_fldtore(dat, EZ);
	mlout_fld(dat, va_N, ha_N, EZ, ez.im()); mlout_fldtoim(dat, EZ);

	mlout_fld(dat, va_N, ha_N, HX, hx.re()); mlout_fldtore(dat, HX);
	mlout_fld(dat, va_N, ha_N, HX, hx.im()); mlout_fldtoim(dat, HX);
	mlout_fld(dat, va_N, ha_N, HY, hy.re()); mlout_fldtore(dat, HY);
	mlout_fld(dat, va_N, ha_N, HY, hy.im()); mlout_fldtoim(dat, HY);
	mlout_fld(dat, va_N, ha_N, HZ, hz.re()); mlout_fldtore(dat, HZ);
	mlout_fld(dat, va_N, ha_N, HZ, hz.im()); mlout_fldtoim(dat, HZ);

	Dmatrix n(va_N, ha_N);
	double dv = (va_end-va_beg)/((double) (va_N-1));
	double dh = (ha_end-ha_beg)/((double) (ha_N-1));
	double vp, hp;

	hp = ha_beg;
	for(int j=0; j<=ha_N-1; ++j)
	{
		vp = va_beg;
		for(int i=0; i<=va_N-1; ++i) { n(i, j) = st.n(vp, hp); vp += dv; }
		hp += dh;
	}
	mlout_n(dat, va_N, ha_N, n);
	mlout_fldviewer(dat, name, va_chr, ha_chr);
	fclose(dat);
	return;
}


/* ------------------------------------------------------------------- 
 * container for field levels at specific points 
 */

/* initialize */
Probe::Probe()
{
	x = 0.0;
	y = 0.0;
	z = 0.0;
	id = '?';
	eps = 1.0;
	ex = CC0;
	ey = CC0;
	ez = CC0;
	hx = CC0;
	hy = CC0;
	hz = CC0;
}
Probe::Probe(double cx, double cy, double cz, char ic)
{
	x = cx;
	y = cy;
	z = cz;
	id = ic;
	eps = 1.0;
	ex = CC0;
	ey = CC0;
	ez = CC0;
	hx = CC0;
	hy = CC0;
	hz = CC0;
}

// field levels
// cp: one of EX, EY, EZ, HX, HY, HZ, SX, SY, SZ
// foa: one of MOD, ORG, SQR, REP, IMP
// id: the field "strength", one of mE, mH, qE, qH, mW 
Complex Probe::field(Fcomp cp) const
{
	Complex f;
	switch(cp)
	{
	case EX: return ex; break;
	case EY: return ey; break;
	case EZ: return ez; break;
	case HX: return hx; break;
	case HY: return hy; break;
	case HZ: return hz; break;
	case SX: 
		f = ey*hz.conj()-ez*hy.conj();
		return Complex(0.5*f.re, 0.0); 
		break;
	case SY: 
		f = ez*hx.conj()-ex*hz.conj();
		return Complex(0.5*f.re, 0.0); 
		break;
	case SZ: 
		f = ex*hy.conj()-ey*hx.conj();
		return Complex(0.5*f.re, 0.0); 
		break;
	default: 
		bundleerror("Probe::field, cp?");
		break;
	
	}
	return CC0;
}
double Probe::field(Fcomp cp, Afo foa) const
{
	return realize(field(cp), foa);
}
double Probe::field(FSid id) const
{
	switch(id)
	{
	case mE:
		return sqrt(ex.sqabs()+ey.sqabs()+ez.sqabs());
		break;
	case mH:
		return sqrt(hx.sqabs()+hy.sqabs()+hz.sqabs());
		break;
	case qE:
		return ex.sqabs()+ey.sqabs()+ez.sqabs();
		break;
	case qH:
		return hx.sqabs()+hy.sqabs()+hz.sqabs();
		break;
	case mW:
		return  0.25*(  MU0*(hx.sqabs()+hy.sqabs()+hz.sqabs())
			+EP0*eps*(ex.sqabs()+ey.sqabs()+ez.sqabs()));
		break;
	default:
		bundleerror("Probe::field, FSid?");
                break;
	}
	return 0.0;
}
   
/* ------------------------------------------------------------------- 
 * an array of probes
 */
   
/* initialize */
ProbeArray::ProbeArray()
          : num(0)
{
	pvec = NULL;
}

/* destroy */
ProbeArray::~ProbeArray()
{
	if(pvec != NULL) delete[] pvec; 
	pvec = NULL;
	num = 0; 
}

/* copy constructor */
ProbeArray::ProbeArray(const ProbeArray& ma)
	  : num(ma.num)
{
	pvec = new Probe[num];
	Probe* ap = ma.pvec;
	Probe* mp =    pvec;
	for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
}
	
/* assignment */
ProbeArray& ProbeArray::operator=(const ProbeArray& ma)
{
	if(this != &ma)
	{
		if(pvec != NULL) delete[] pvec;
		num = ma.num;
		pvec = new Probe[num];
		Probe *ap = ma.pvec;
		Probe *mp =    pvec;
		for(int i=0; i<=num-1; ++i) *mp++ = *ap++;
	}
	return *this;
}

/* delete all Probe entries */
void ProbeArray::clear()
{
	if(pvec != NULL) delete[] pvec;
	pvec = NULL;
	num = 0; 
}

/* subscripting */
Probe& ProbeArray::operator() (int i)  
{ 
	if(i<0 || i>=num) bundleerror("ProbeArray: () out of range");
	return pvec[i]; 
}
Probe ProbeArray::operator() (int i) const 
{ 
	if(i<0 || i>=num) bundleerror("ProbeArray: () out of range");
	return pvec[i]; 
}

/* add a probe */
void ProbeArray::add(Probe m)
{
	Probe *avec;
	avec = new Probe[num+1];
	Probe* ap = avec;
	Probe* mp = pvec;
	for(int i=0; i<=num-1; ++i) *ap++ = *mp++;
	*ap = m;
	if(pvec != NULL) delete[] pvec;
	pvec = avec;
	++num;
	return;
}